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Adaptive Neural Controls (ANC) Executive Summary 

Despite recent advances in efficiency, current methodologies for space structure control design 
still engage significant human resources for engineering development and maintenance. For 
example, since fixed-gain space system controllers must be updated periodically to adjust to in¬ 
mission changes in system dynamics, this implies burdensome ground support activities. 

The Adaptive Neural Controls (ANC) program is part of an effort to develop neural network 
based controllers capable of self-optimization, on-line adaptation and autonomous fault detection 
and control recovery. This development also supports the Nation’s long-term space exploration 
objectives for which autonomous spacecraft involving self-reliant control systems are a necessity. 

The ANC Program was funded in two phases. The first, Basic phase focussed on the 
development of efficient and completely autonomous neural network feedforward control for the 
case of broadband disturbances. Algorithms were developed that work with no prior modeling 
information about the system to be controlled and adapt to changing conditions, while 
minimizing or eliminating the introduction of extraneous training signals. The algorithms were 
demonstrated experimentally on an optical structure testbed at Harris. 

The second, Amended Option 1 phase of the program demonstrated a more complex neural 
controller on the ASTREX test facility at the Air Force Phillips Laboratory capable of fault- 
tolerant adaptive control of multiple sensors and actuators. This system used 6 actuation 
channels of the existing ACESA struts on the ASTREX structure to simultaneously cancel three 
independent tonal disturbances in the 10-15 Hz band, measured at non-collocated sensors on the 
secondary tower of the structure. The system demonstrated impressive fault-recovery 
performance, maintaining good cancellation performance with successive actuators disabled, 
down to a minimum of two out of six. Cancellation of individual tones was between 20 and 50 
dB, with over 28 dB attenuation realized RMS. The algorithm required very low computational 
throughput, operating at a sample rate of 1/20 Hz. 

The results of the ANC program show that adaptive cancellation systems can reduce vibrations in 
precision structures without prior modeling data and can adapt successfully to failures in 
actuators or sensors, optimally reconfiguring themselves without human intervention. These 
capabilities should significantly reduce the expense of designing and maintaining vibration 
control systems for spacecraft. 

The next step in the development of ANC technology should be to demonstrate broadband 
adaptive neural cancellation using a small flight-like processor such as the Modular Control 
Patch. This would give potential users confidence that the technology could be integrated into 
their existing or planned systems with little overhead in terms of cost, volume, weight, or power. 
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1. Introduction 

With the advent of Adaptive Neural Control (ANC), we are on the threshold of revolutionary 
changes in control design and implementation. These changes will impact not only structural 
vibration control but a wide variety of other controls tasks as well. Just in the area of vibration 
control, the advances made recently in ANC translate into greater cost-effectiveness for control 
design, implementation and verification and a corresponding expansion of DoD and commercial 
applications. 

By way of introduction, we first trace the historical trends in active structural controls that have 
stimulated and provided the basis for ANC development at Harris. We devote careful effort to 
this discussion in order to convey our vision of how ANC fits in with the on-going evolution of 
structural control technology. The historical record suggests that any technology undergoes 
several distinct phases in its development: 1) an embryonic phase in which fundamental data is 
gathered and the essential difficulties are identified and grappled with, (2) a basic competence 
stage in which full-scale capabilities are demonstrated but the technology is still costly and 
cumbersome in its implementation and not necessarily effective in its results and, finally (3) the 
incisive stage wherein not only can the technology be applied routinely with the desired success 
but the processes of its application are quick and inexpensive. The third, incisive stage must be 
reached in order for the technology to find use in a broad range of industrial and commercial 
applications. 

The area of active structural control design for aerospace and commercial systems already display 
at least the first two development phases noted above. Although, at present, there are groups still 
working within nearly every point of the development sequence, the state of the art is now 
entering the third phase. 

The first, embryonic phase of structural control design began in the mid 70’s and extended into 
the early 80’s. This early stage involved much theoretical effort in its attempts to apply "modem 
control theory" (Kalman filtering and optimal feedback control theory developed in the 60’s) to 
vibration control problems. During this period also, many sub-scale, highly idealized 
experiments were performed to test the basic concepts and identify the significant roadblocks to 
effective performance. As a result, several difficulties with the earlier control theories were 
identified and efforts were made to resolve them. For example, two of the most severe 
challenges were the design of (1) simplified controllers (low order dynamic controllers for high 
order systems) and (2) robust controllers that maintain both stability and acceptable performance 
in the face of system modelling errors and in-mission changes in physical characteristics. In 
response to these challenges, a variety of reduced-order and robust control design methods were 
developed and the feasibility of at least a subset were demonstrated via small scale laboratory 
experiments. Although this first development phase produced most of the fundamental tools and 
illustrated their feasibility, convincing evidence was still lacking for the claimed performance 
benefits of the technology when applied to full scale structural control problems. 

The second, basic competence, phase of development involved the competitive sorting out of a 
plethora of various design approaches, using test experiences on larger scale, nonidealized 
apparatus. A really positive stimulus was provided by the emergence of "Guest Investigator" 
efforts which provided several investigative groups the opportunity to conduct research on 
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traceable testbeds that would otherwise be beyond the technical and financial means of individual 
researchers. An example is the NASA Controls/Structures Integration (CSI) Guest Investigator 
(GI) Program, that recently completed its second phase. Both guest investigator programs and 
other large scale experimental efforts involved independently refereed testbeds—i.e., the 
hardware is provided and hosted and the control algorithms are actually implemented (and test 
data taken) by an independent Government appointed body. This arrangement provided a 
completely independent mechanism for unbiased, critical review of competing approaches and 
we restrict discussion here to independently referred test results only. Also, of course, our 
discussion includes only those efforts in which Hams researchers participated. This experience 
however, amply illustrates how the basic competence phase of technology development was 
successfully negotiated and how the way was opened to the third stage of development. 

Harris researchers have demonstrated vibration control design on six independently refereed 
testbeds. In Figure 1-1 we show, for each testbed, the values of: 

Number of Performance Degrees of Freedom 

= number of generalized degrees of 
freedom that must be independently 
suppressed to address performance 
requirements 

Number of Modes in Disturbance Band 

= number of vibration modes in the 

disturbance band that must be attenuated 
by at least an order-of-magnitude in order 
to meet performance requirements 

These two variables give a rough indication of the difficulty of the control problem. The ACES, 
Mini-MAST, CASES, CEM Phase 2, and CEM Phase 3 test facilities were addressed on the 
NASA CSI GI program in the order given. Harris was the only investigator group to participate 
on both phases of this program. The performance results achieved are documented in [1]. Most 
of the testbed experiments involved only control algorithm design with the exception of CEM 
Phase 3 and ACTS which also involved test demonstrations of Harris' inertial proof-mass 
actuators and active isolation intrastructural actuators. 

The above testing experiences using large scale facilities under nonidealized conditions, helped 
us (and other groups) progress beyond the embryonic stage of development in two important 
ways. First, the early experimentation on ACES and Mini-MAST showed that the development 
of advanced design tools (such as Matlab toolboxes for Maximum Entropy/Optimal Control, 
p-Synthesis, or Quantitative Feedback Design) is not enough for efficient (or even successful) 
implementation of working control systems in practice. What is needed, in addition, is a 
practical methodology for using advanced design tools that meshes algorithm design, system 
modeling and subsystem and component tests, into a realistic strategy for implementation. The 
development of such a methodology was key to our success. The second area of progress was the 
streamlining of the whole design/test/validation process to achieve validated control results with 
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less time and expenditure. For example, the cumbersome, time consuming process of building 
and refining finite-element models was replaced by the use of automated system identification 
methods for extracting dynamic models directly from test data. This approach fits in with on- 
orbit control refinement for space systems or with any kind of "retrofit" situation in which 
vibration control is desired for existing structural hardware. A second area of streamlining was 
the introduction of faster design tools that permitted design turnaround on-site. 



No. of Modes in the Disturbance Band 


Figure 1-1: Six Independently Refereed Vibration Control Experiments have been 
Completed by Harris Researchers 

The development of a systematic design/test methodology and its subsequent streamlining 
significantly improved efficiency over the course of the GI program. For example, the ACES 
facility at NASA/MSFC (the first testbed addressed) required a man year and three on-site test 
sessions (2 days each): An initial session for open-loop testing, a second to test an initial, robust 
controller (while collecting closed-loop ID data) and a third session to test the integrated, 
advanced controller. In contrast, a few years later, the CEM Phase 2 testbed, required two 
months of design effort and a single on-site test session. This result is obtained despite the fact 
that the CEM Phase 2 and ACES facilities were of comparable difficulty (see Figure 1-1). 
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The above efficiencies were essential to successfully address the most challenging structural 
control problem yet posed: the ACTS program. 

On this program, the structural testbed (Figure 1-2) is provided by an independent host 
contractor. The testbed is very complex, traceable and program-directed. Disturbances at the 
base have flat PSDs from 10 to 200 Hz. Stringent performance requirements are levied on both 
line-of-sight (LOS) and aberration errors, such as defocus and primary reflector surface 
wavefront errors. Because of the very broadband disturbances, tens of modes contribute to LOS 
error and more than a hundred modes contribute to wavefront error. In all, the response of over 
150 modes must be attenuated by more than 20 to 40 dB in order to meet performance 
specifications. The Harris effort involves not only control algorithm design but also the 
fabrication, integration and test of two distinct hardware complements: 24 Advanced LPACTs 
and 6 Active Isolator Fittings (AIFs). In all, the control design problem entails 30 inputs and 30 
outputs. Our mission was to quantify, as comprehensively as possible, the complexity and cost 
versus performance tradeoffs inherent in the application of active structural control to this class 
of systems. 



Test data shows over 40 dB attenuation of LOS related modes. Very significant results were also 
achieved for aberration errors. Figure 1-3 shows open-versus closed-loop PSD test data results 
for defocus. Nearly 50 dB attenuation is achieved for the dominant modal peak and over 30 dB 
rms is observed over the 10 to 200 Hz frequency band. Finally, Figure 1 -4 shows open-versus 
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closed-loop PSD data for wavefront error. Evident in this plot is the large number of modes 
contributing to wavefront error in the open-loop. Dominant modal peaks are reduced by nearly 
40 dB and over 20 dB rms attenuation is obtained over the entire disturbance band. 

In all of the test programs in Figure 1-1, the fixed-gain control designs invariably worked on the 
first trial and performance specifications for each testbed were invariably met or exceeded. 
Because of the effectiveness and predictability of the results, the independently refereed testing 
experiences did much to make the technology an area of professional competence backed by 
reliable, effective tools —and this is the worthy goal of the second, "basic competence” phase of 
technology development. Our efforts to streamline the design and test process not only reduced 
cost and engineering design time but also led us to the threshold of the third, "incisive" phase of 
development. Our results suggested the idea that if one can systematize and streamline a human- 
operated design/test process (such that modeling, control design and test verification are all 
accomplished in a single on-site session), then one ought to be able to automate the entire 



Figure 1-3: ACTS Defocus Errors are Reduced by More Than 30 dB rms Over the 
10-200 Hz Disturbance Band and Peak Amplitudes are Attenuated by Nearly 50 dB 
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Figure 1-4: ACTS Wavefront Errors are Reduced by More than 20 dB rms Over the 10- 
200 Hz Disturbance Band and Peak Amplitudes are Attenuated by More than 40 dB 

process. This is the first motivation for ANC. The second motivation arose because the 
invention of a new neural network architecture for addressing dynamic systems by Dr. D.C. 
Hyland in 1989 provided an algorithmic means for realizing the desired automation. Finally, 
over the intervening years, advances in Digital Signal Processor (DSP) technology provided ever 
faster and cheaper hardware capable of implementing the ANC algorithms in real time. 

As mentioned, the incisive stage of technology development is attained when not only is the 
technology proved capable of reliable results but the implementation of the technology is rapid 
and relatively cheap. Total automation of the control design/test process is needed because the 
current methodology still engages significant human resources. Since space system control 
designs involving fixed-gain controllers must be updated periodically to adjust to in-mission 
changes in system dynamics, this implies burdensome ground support activities. Besides 
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reducing engineering manpower requirements, such advances in automation support the Nation’s 
long-term space exploration objectives for which autonomous spacecraft involving a self-reliant 
control systems are a necessity. Such robot explorers would have to independently update 
control laws, detect faults and reconfigure control systems. Finally, while there are numerous 
commercial areas that would use or benefit from active vibration control, this will not occur until 
dramatic decreases in engineering development cost and dramatic improvements in system self- 
reliance are achieved. 


7 









ANC Final Report 


2i The Adaptive Neural Control Architecture—Overview 
Description 

For the past four years, motivated by the technology needs described above, Harris has been 
developing a new neural architecture to implement on-line systems identification and adaptive 
control for systems. Applications were originally targeted to vibration suppression in precision 
space structures but have since been expanded to noise and vibration control and motion control 
for both DoD and commercial systems. Our studies began in 1989 with the discovery by Dr. 

D.C. Hyland of a new neural architecture for identification and control. The new architecture, 
termed the Adaptive Neural Control (ANC) architecture, subsequently led to a sequence of 
successful demonstrations and new development efforts. 

Much of the previous work in adaptive control via neural networks (see [2] for an excellent 
review) concentrated on highly nonlinear but low dimensional systems. In contrast, the ANC 
architecture concentrates on neural schemes particularly geared to problems involving high order 
systems exhibiting very broadband dynamics. As indicated in Figure 2-1 ANC combines tapped 
delay lines with "static" neurons (each neuron is a two-way device incorporating a back 
propagation path) to perform on-line system identification and adaptive control. The system 
adapts in the presence of unknown persistent plant disturbances and instrumentation noise and 
requires no detailed prior modeling information. 

There are several key features of this architecture that have made it particularly attractive. First, 
although the architecture can be visualized as a neural network, the control scheme is 
fundamentally a massively parallel, decentralized adaptive control algorithm that need not be 
implemented literally as a collection of artificial neurons. Secondly, these "neural" algorithms 
feature learning capability that is distributed down to the smallest computational unit. 
Decentralization (distributed learning) imparts the ability to autonomously recover from 
hardware failures—including damage to the neural processor itself. A third key feature is that 
the basic neural building blocks are hierarchically organized into a set of standardized modules. 
Analogous to a "Lego set," modules can be combined to build an enormous variety of systems 
and permits complex systems to be built up from simpler components in a transparent way. 

Finally, modularity and parallelism yield implementation flexibility. Specialized hardware is not 
required for implementation of the Harris ANC architecture. The entire identification or control 
algorithm can be distributed among several parallel processors, and hardware suitable for this 
purpose is currently available and is being used for engineering development. This means that 
we can progress in orderly fashion from the use of existing ICs to (ultimately) dedicated neural 
ICs, thereby building our capabilities gradually and systematically. 

While details of the ANC architecture are given in recent papers and reports [3]-[5], we briefly 
review the basic features here. The hierarchy of modular structures is shown in Figure 2-2. This 
hierarchy starts, at the lowest level, with tapped delay lines and neurons with intrinsic back 
propagation. These are the same "static" neurons that would be utilized for such applications as 
pattern classification and nonlinear mapping. The key to applying such neurons to dynamic system 
identification is to organize them into larger building blocks, the dynamic ganglia. A ganglion is 
an array of neurons designed to establish temporal ordering within the network so as to process 
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time histories of network signals. Ganglia are interconnected by bundles of synapses, which are 
sometimes constrained so that their weights form Toeplitz matrices. 



Figure 2-1: The Adaptive Neural Controller Executes Simultaneous System Identification 
and Adaptively Optimized Control 

The next level in the hierarchy combines ganglia and their connecting synaptic bundles to form 
replicator units. The basic job of a replicator unit is to duplicate the output of a previously 
unknown sampled - data dynamic system when both replicator and system are stimulated by the 
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same training input. Thus the replicator is the basic module for system identification. Several 
types of replicator have been developed, each corresponding to a particular model form. The 
work in identification methods using systems observer Markov parameters [6], [7], led to 
discovery of a new model form for dynamic systems - the ARMarkov model, so-called 



Figure 2-2: Hierarchy of Modular Nebular Structures Progressing from Basic 
Constituents to Higher-Level Modules 

because it combines features of impulse response (Markov parameters) with ARM A (Auto 
Regressive Moving Average) models. 

Furthermore, most of the replicator types have both time and frequency-domain (or mixed time 
and frequency-domain) versions. Generally, one obtains a frequency-domain version of a 
particular time-domain replicator by replacing the time series inputs and the time series error 
signals by the outputs of Finite Fourier Transform (FFT) filters having these signals as inputs. 
This operation is equivalent to inserting a multiplication by a unitary matrix within the replicator 
structure. However, this superficially trivial change allows the replicator to address frequency 
weighted output matching and direct frequency response identification. 

Numerous analytical examples have been produced to demonstrate structural identification or 
control using neural replicators. Some of these involve the use of simulated input/output data 
and others use actual test data. Also, a number of laboratory experiments have been performed. 
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In some cases a MATLAB simulation was used to implement the neural algorithm, while in other 
cases the algorithm was implemented in real time using a DSP card. For example: 

1. Figure 2-3 shows a simple beam experiment that produced excellent convergence of the 
adaptive model to the actual structural plant in 100 seconds. 

2. The neural network system identification capability was also demonstrated on the Harris 
Multi-Hex Prototype Experiment (MHPE), which is a four meter Cassagrain test 
structure. Figure 2-4 shows the network converged to the MHPE plant in 125 seconds. 

3. Using a Digital Signal Processor (DSP), an Internal Research and Development (IR&D) 
experiment in active acoustical noise cancellation was completed in which over 20 dB 
broadband attenuation was achieved. 

Many of these examples involve multiple inputs and outputs and nearly all involve fairly 
complex structures with many modes in the frequency band of interest. Also the laboratory 
experiments tested the algorithm under such real-world complications as sensor noise and 
ambient steady-state disturbances. Summarizing this experience, we can say that reasonably 
complex multi-mode systems can be identified with excellent accuracy with convergence times 
ranging from a few minutes to a few seconds (depending on numerous factors, such as system 
sample rate, frequency band of interest, etc.). 

Returning now to the hierarchy shown in Figure 2-2, several replicator units are combined in 
order to form the Adaptive Neural Control (ANC) system. An ANC performs on-line, 
simultaneous system identification and adaptively optimized control. The most basic ANC 
architecture for simultaneously replicating an unknown plant and adapting an output feedback 
controller so as to match the closed-loop input/output characteristics with a prescribed reference 
system has two parts: (1) the closed-loop modeller and (2) the control adaptor. The closed-loop 
modeller uses training signals and the plant sensor output to adapt the weights so that the closed- 
loop is replicated. After convergence, the modeller output matches the closed-loop system - in 
effect the modeller identifies the plant within the closed loop. 

In the control adaptor, there is an internal model of the plant, copied from the plant modeller. 

Thus the control adaptor can, in effect, back-propagate error through the plant model to the 
controller output location. With its internal model of the plant, the adaptor uses the training 
signal, its own output and that of the reference system to adjust its weights so that the reference 
system is replicated. 

The above form of ANC addresses the feedback type of control architecture. It is useful here to 
note the two quintessential architectures for adaptive control: feedforward cancellation and 
feedback control. Figure 2-5 illustrates these two basic disturbance suppression strategies. In 
feedforward cancellation, the system motion measurements (the error sensor outputs) are used 
only to adjust the feedforward gain, F. Once F converges, the actuators are driven by the 
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Mini-Mast Testbed Facility Ajjg,Performance Goals: Reduce firs! betiding 

mode response by 
> 20 dB 



Figure 2-3: Beam Experiment Shows Excellent Convergence Results After 100 Seconds 
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Multi-Hex Prototype Experiment (MHPE) 
Addresses Vibration In Large Optics 



4M Diameter Cassagrain 
Configuration 

Vibration Control System 
Uses Nine LPACT Sensor/ 
Actuator Units 


MHPE System ID Results 

— Measured FRF 

- Neural Replicator 

Frequency Response 



Figure 2-4: Neural Network Closed-Loop System Identification Experiment was 
Performed on the Dynamically Complex Multi-Hex Prototype Experiment (MHPE) 
Testbed 
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(♦denotes convolution) 

Once F converges, m directly drives u 


Adaptive 

Feedforward 

Cancellation 



Adaptive 

Feedback 

Control 


* Comparable or 
superior performance 
capabilities 

* Better transient 
response 

* Can simplify 
control hardware 

* Can handle disturbances 
that are inaccessible 

to measurement 


Figure 2-5. Feedforward vs. Feedback Architecture 
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disturbance measurements (assuming they are available) or a synchronization signal or timing 
pulse (if the disturbances are periodic). This control architecture is workable only if the 
fundamental disturbance sources are accessible to measurement. In contrast, adaptive feedback 
uses the motion measurements not only to adjust the feedback gain B but also to drive the 
actuators. This is the suitable architecture for systems wherein disturbances are not accessible to 
direct measurement. Adaptive feedback is a much more challenging problem than feedforward 
and its successful treatment by ANC is a significant advance. Moreover, adaptive neural controls 
have been elaborated into a selection of forms capable of handling feedforward or feedback or a 
combination of the two. This last feature is unique. 

In addition, the above time-domain versions of ANC also have frequency-domain counterparts 
(combining several frequency-domain replicators). In the frequency-domain ANC that combines 
feedforward and feedback, the two forms of control have synergistic benefits. First feedback 
control is able to attenuate unmeasurable disturbances and can greatly improve transient 
response, thereby permitting faster feedforward convergence. Also, closed-loop feedback 
smoothes out the system transfer functions. This allows the feedforward control to track much 
more rapid system changes—e.g. disturbance frequency variations due to motor 
acceleration/deceleration. Alternately, smoother transfer function variations permit a simpler and 
cheaper feedforward system to be used. 

The earliest ANC efforts concentrated on feedback control and many analytical examples were 
investigated. 

Figure 2-6 shows an early example of ANC operation. Using a simulation model of the Mini- 
MAST facility, an ANC simultaneously performed system identification and control 
optimization. In this example, the ANC was required to achieve more than an order of 
magnitude closed-loop attenuation of the first bending mode pair of Mini-MAST leaving higher 
frequency modes unaltered. This basic control objective was obtained within 7.5 sec. of 
adaptation and exact agreement with desired closed-loop response was attained after four 
minutes of adaptation. 

In its detailed operation, the ANC carries out a sequence of steps analogous to the modeling ID 
and design refinement steps carried out by human designers within a streamlined design and test 
methodology. However, the ANC carries out these steps tremendously faster and without direct 
human supervision. 

The above example helps to illustrate the potentially enormous savings in time and effort for 
development of initial space structure control design. Note that, by typical performance 
standards set by the human G.I.s on the Mini-MAST testbed during the CSI GI Program Phase 1, 
the control design obtained by ANC (Figure 2-6) is quite respectable. However, rather than 
requiring an elaborate design model together with a man year of effort with several hours of on¬ 
site testing, ANC obtains its results without prior information on the testbed and within less than 
five minutes of unsupervised operation! 
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Neural Network Experimental Setup 



• Beam is cul from .125 inch (hick tluminun sheet 

• Beam is cantilevered horizontally 

• LPACT acts in Z-aais (normal to paper) 


• Training signal is logarithmic sine sweep from 5-40 Hz. 

• 180 seconds of data stored 

• A/D‘s docked at 80 llz sample rate 



Frequency Response Measurement 
Neural Replicator Frequency Response. 


Figure 2-6: Simulation Results for the LaRC Mini-MAST Testbed Demonstrate 
Simultaneous System Identification and Control Optimization 

The above examples stimulated interest in demonstrating the ANC in the laboratory. Recently, 
feedback control was demonstrated by NASA/LaRC in collaboration with Harris using the 
Controls Evolutionary Model (CEM) testbed at NASA/LaRC in July 1993. A simplified version 
of the general (time-domain) ANC model reference adaptive control scheme was used, in which 
the on-line controller was constrained to be an optimal one-step-ahead controller for the 
identified plant. The entire algorithm was implemented by programming the lab facility 
computer. Figure 2-7 summaries the test results. Broadband inputs were used as training stimuli 
to execute identification and control adaptation. With approximately 10 seconds of training, the 
system converged to an optimal deadbeat controller. 

The above results indicate at least the feasibility of ANC and provided valuable experimental 
experience with which to avoid potential pitfalls. On the other hand, the experimental tests made 
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little use of the highly modular, parallel character of the ANC architecture—the algorithm was 
programmed as a monolithic entity in the lab facility computer, which had capabilities far in 
excess of that needed. Sizing studies of ANC algorithms have shown that, judging from 
throughput and I/O bandwidth requirements, the most general ANC can be implemented in a 
parallel system of several DSP cards, using existing hardware. This configuration is more 
suitable for practical DoD applications and is a necessity for commercial or dual-use 
applications. 

The development and testing of such DSP-based implementations is the subject of the next 
chapter. 
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3. Demonstration of Frequency-Domain ANC Algorithms for 
Feedforward Control 

Following the fundamental theoretical and experimental investigations described above (in which 
the basic feasibility of Adaptive Neural Control was established), development efforts focused on 
implementation using practical, productizable hardware. It was concluded that a reasonable, 
near-term hardware implementation suitable for spacecraft and dual use applications involves the 
use of existing DSP hardware in a multi-processor architecture. Consequently, hardware of this 
type was used to implement ANC and intensive experimental tests were conducted on a variety 
of testbeds. 

Much of the DSP experimentation was carried out under a Corporate funded project for 
commercialization of the ANC technology. Thus, this effort concentrated on developing ANC 
algorithms for DSP hardware that could cost-effectively support commercial noise suppression 
and vibration control applications. A survey of these applications reveals that the most useful 
distinction from the point of view of system complexity is the spectral distribution of 
disturbances. As Figure 3-1 illustrates, there is a hierarchy of disturbance types ranging from 
harmonic or "single-tone" disturbances, which are the simplest to address, to "broadband" or 
continuous spectrum disturbances which constitute the most general case and are the most 
complex. 

The harmonic disturbance case refers to those disturbances that are dominated by a single 
sinusoid, where the amplitude and frequency may vary over time (at rates somewhat below the 
dominant frequency). It appears that disturbance sources for the majority of commercial 
applications fall into this class because the objectionable noise or vibration is directly caused by a 
rotary motor or engine. Examples include fan noise, automobile engine noise, aircraft cabin 
noise, etc. In many of these cases, the disturbance source—i.e. the engine, motor, rotor or CMG, 
etc., is directly accessible to measurement, either through a synchronization or timing pulse pick- 
off directly from the motor or by use of motion sensors located on or very near the source. 
Consideration of the overall system dynamics shows that the knowledge that the disturbance is 
harmonic combined with an in-situ timing pulse measurement is tantamount to a direct 
measurement of the disturbance forces. Consequently, the feedforward control architecture is 
suitable for suppression of motion due to these disturbances. 

At the next level of complexity, multi-tone disturbances consist of two or more harmonics, where 
again the frequencies and amplitudes of each component may be subject to time variation. This 
disturbance type is significant in most applications, e.g., nonlinearities in engines and motor 
mounts produced higher harmonics and sub harmonics in addition to the fundamental tone 
produced by the rotary source. In many cases, the multi-tone disturbance is truly periodic —i.e. 
all constituent tones are multiples of the fundamental frequency. There are other cases where the 
constituent tones are incommensurate —i.e. although multi-tone as defined above, the disturbance 
as a whole is periodic. This occurs for many systems of DoD interest. For example, a spacecraft 
bus many contain several Control Moment Gyros (CMGs) for attitude control Each CMG, due 
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• Most applications involve both multi-tone and broadband to some degree 

• A complete capability must address all 3 types of disturbances 

• The three disturbance types form a hierarchy of increasing difficulty 



Non-negligible component in many applications. 
Significant component in a subset of applications. 
Example: Dishwashing, washing machine noise 
has a big component due to fluid sloshing (which is 
broadband) Also: Droadband aerodynamic noise in 
aircraft cabins 


Significant (but not dominant) feature in most 
applications -e.g. nonlinearities in engines and 
motor mounts produce higher harmonics besides 
a fundamental harmonics 


MAIN component inthe largest number of 
applications - e.g. autombile engine noise, aircraft 
cabin noise, fan noise, etc. Basically any instance 
of engine/motor driven disturbances 


Figure 3-1: Applications of Adaptive Control Involve a Hierarchy of Disturbance Types 
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to static and dynamic imbalances produces a primarily harmonic disturbance at its rotor speed. 
Although nominally identical, each CMG has a slightly different rotor speed. Thus the 
vibrational disturbances imparted to a sensitive payload mounted on the bus consists of several 
harmonics with closely spaced frequencies. The periodic case with closely spaced frequencies is 
particularly difficult to treat and has been the focus of our ANC developments for multi-tone 
disturbances. Obviously, as in harmonic disturbances the fundamental rotary disturbance sources 
are also accessible to direct measurement—i.e. one can utilize synch signals taken off each rotary 
source, etc. In addition, however, in several DoD applications, one may be constrained to use 
ambient motion measurements that are not proximate to the disturbance sources. The successful 
treatment of such cases has been an objective of the ANC development. Whether the disturbance 
measurements are obtained directly or indirectly from ambient motion sensors, the feedforward 
control architecture is the most suitable approach for multi-tone disturbances. 

Finally, at the highest level of complexity, broadband (or "continuous spectrum") disturbances 
are defined as having continuously distributed power spectral densities (PSDs), without discrete 
tones. This is a nonnegligible component in many applications; e.g. fluid sloshing noise in dish 
washers and washing machines. This is also a nearly dominant component in a subset of 
applications; e.g. broadband noise in aircraft cabins produced by aerodynamic turbulence. In 
many cases, these disturbances are not accessible to direct measurement so that th tfeedback 
control architecture must be used. In other cases, even though the disturbances are amenable to 
(at least indirect) measurement, the broad spectral distribution of the disturbances make a 
combined feedforward and feedback control architecture the most efficient, robust approach. The 
ability to implement adaptive feedback and combined feedforward, feedback control 
architectures is a unique strength of ANC technology. An example of broadband disturbances 
that are accessible to measurement is the case of vibration in a large optical structures due to 
rapid repointing maneuvers—a situation addressed by the Phillips Lab ASTREX testbed. Here 
the ultimate disturbance sources are the command signals to the slewing actuators—these are 
generated by the onboard computer and are obviously directly accessible to the vibration 
controller for combined feedforward/feedback control. 

On the Corporate funded commercialization projects, we implemented a variety of ANC 
algorithms for harmonic and multi-tone disturbances using a PC interfaced DT3808 DSP card (TI 
TMS320C40 chip) and developed efficient, rapidly convergent adaptive controls via real-time 
testing on three entirely different laboratory testbeds. These testbeds are (1) an acoustic duct 
apparatus, (2) a concentrically mounted proof-mass actuator test rig and (3) the MHPE optical 
structure vibration experiment apparatus. We first describe the testbeds used as follows: 

The Acoustic Duct: 

Figure 3-2 shows an engineering drawing of the apparatus and Figure 3-3 shows a simplified 
block diagram of the testbed. The duct is a 20 ft. long tube with 8 in. inner diameter. It is 
mounted at one end with a disturbance speaker and the actuator speaker can be placed at two 
alternate locations. The objective is to suppress noise received at the end of the tube as 
monitored by an acoustic microphone. 
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Figure 3-2: The Acoustic Duct Apparatus was Used to Investigate Adaptive Control of One 
Dimensional Acoustic Systems 
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alternative site 



Figure 3-3: Schematic Diagram of the Adaptive Control System for the Acoustic Duct 
Apparatus 
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Concentrically Mounted Actuator Test Rig 

The concentrically mounted proof-mass actuator test rig is shown in Figure 3-4. A base plate is 
connected by a tubular structural member to a star-shaped test body (at the top of Figure 3-4) 
having nontrivial dynamic characteristics. The system is disturbed by up to three shakers 
(labeled dist #1 etc. in Figure 3-4) mounted on the base plate. A proof-mass actuator, having a 
toroidal configuration is clamped around the lateral surface of the tube in order to provide control 
forces. The objective is to suppress vibration in the star shaped body at the top of the test rig, as 
monitored by an accelerometer at the top of the tube. 

MHPE Apparatus 

The optical structure testbed is the Harris Multi-Hex Prototype Experiment (MHPE) apparatus, 
constructed under a Harris IR&D program in 1988. Figure 3-5 shows a simplified block diagram 
of the MHPE and the test setup. The MHPE has the form of a Cassagrain telescope and consists 
of a base plate (supported by air bags) connected by a six member truss to the primary mirror 
support structure, holding a tripod tower which supports the secondary mirror assembly. 
Disturbances are injected using a proof-mass actuator affixed to the primary mirror support 
structure. Three additional proof-mass actuators with colocated accelerometers are mounted on 
the secondary mirror platform to provide actuation and motion sensing. The objective is to 
suppress vibrational accelerations at the secondary mirror platform that are excited by the forces 
applied to the primary support structure. 

For feedforward cancellation of single-tone disturbances in all of the above three testbeds, we 
specialized the feedforward, frequency-domain version of the ANC algorithm. In contrast with 
basic LMS algorithms, the ANC system for this case is able to simultaneously identify all needed 
transfer functions and simultaneously adapt the actuator inputs without interrupting normal 
operation or injecting an extraneous, broadband test signal (dither). Moreover, the ANC 
algorithm is fast: with no previous identification, completed cancellation is achieved in 3 to 5 
iterations (10-15 seconds in these experiments); once transfer coefficients are identified, control 
adjustment can be accomplished in one step. 

The single-tone ANC algorithm was experimentally demonstrated on all three testbeds. These 
live demonstrations are documented in the Harris Corporation video "Adaptive Noise and 
Vibration Cancellation Demo". Starting from a clean slate (no prior identification data) the 
system converges in a fraction of a second. The video shows how ANC can track smooth 
changes in system dynamics, and can quickly recover after dramatic, sudden system changes. 
Moreover, the same algorithm, without modification is shown to work on many different types of 
systems. Over 20 dB noise or vibration reduction is achieved in all cases. 
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Figure 3-4: Force Cancellation System Demonstration 
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Diagnostic 



Figure 3-5: System Block Diagram for Multi Hex Prototype Experiment (MHPE) 
Demonstration of Harmonic Vibration Cancellation 
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Subsequeptly, the ANC single-tone algorithm was extended to handle multiple tones. The ANC 
system for this case combines a fast neural demodulator unit with a array of single tone canellors. 
Multi-tone cancellation is routinely demonstrated in the laboratory. For example. Figure 3-6 
shows results pertaining to the concentrically mounted actuator test rig for two relatively closely 
spaced harmonics. Over 30 dB and approximately 40 dB attenuation is achieved on the two 
tones, respectively. Starting with no prior transfer coefficient information, simultaneous 
identification and control are achieved in approximately one second. For the same testbed. 

Figure 3-7 illustrates results for three closely-spaced tones. The algorithm again converges in 
approximately one second, achieving 20 dB attenuation for all tones. 

The above demonstrations revealed some unique features in our ANC, frequency-domain 
algorithms. As noted, one unique feature is the simultaneous control adaptation and system 
identification without test signals and without the need for look-up tables. Although this was not 
done on the Corporate-supported projects, the potential of ANC to combine feedforward with 
feedback control is also unique. 

With the above background, we proceeded on the ANC program to investigate feedforward 
control for continuous spectrum disturbances—considering both time-domain and frequency- 
domain versions of ANC. 
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Figure 3-6: Multi-Tone Neural Algorithm Results for the CMLPACT Test Rig 
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4. ANC Program Developments 


4.1 Introduction 

The Basic segment of the ANC program is essentially a preparation for experimentation on the 
ASTREX testbed (or a suitable alternative) with the objective being to suppress vibrations 
excited by repointing maneuvers. These disturbances are typically continuous spectrum, not 
tonal. As noted for the repointing situation in the previous chapter the ultimate disturbance 
sources are the command signals to the slewing actuators and these are directly accessible to the 
vibration controller. Accordingly a feedforward control architecture is appropriate. Hence on the 
ANC Basic program we concentrated on feedforward control of continuous spectrum 
disturbances. We studied various issues associated with the use of DSP hardware using the 
MHPE testbed. The objective of this phase is to successfully demonstrate adaptive control of the 
MHPE testbed with an ANC algorithm that can be readily scaled up to handle larger scale test 
facilities. 

So that we may give a more detailed description, Figure 4.1-1 again shows the test facility used 
to-date. The MHPE without its outboard panels is mounted on air bags in the inflated condition. 
A Linear Precision Actuator (LPACT), proof-mass actuator is mounted on the center panel to 
provide broadband disturbances emulating the vibrational disturbances that might be excited by 
equipment or rapid slewing. The control objective is to suppress lateral motion of the secondary 
mirror platform since this is directly proportional to line-of-sight (LOS) errors. Accordingly, 
only the LPACTs on the secondary mirror platform are used for control and tower bending 
motion is sensed by the colocated accelerometers. These actuators and sensors interface with a 
PC based DT 3808 DSP card, as does a spectrum analyzer-generated disturbance source. 

Using this test set-up, two main issues were explored. Although, in the tonal cancellation work 
described in the proceeding chapter, a frequency-domain version of ANC was the obvious 
choice, both time-domain and frequency-domain forms can be used for continuous spectrum 
disturbances. Thus our first objective was to evaluate both time and frequency-domain versions 
of ANC algorithms. Secondly, while throughput, I/O bandwidth and memory capacity 
requirements have been carefully studied for general ANC algorithms, the design limitations 
imposed by instrumentation noise and machine accuracy; i.e. word length, needed further study. 
Thus, a second objective was to shed more light on these issues by means of extensive 
experimentation with the DSP hardware. 

Our approach in conducting these studies was to closely enmesh theoretical deductions with 
numerous and frequent experimental trials. Table 4.1 summarizes in chronological order the 
thirteen major experimental tests that were performed using both time and frequency-domain 
versions of ANC, before the recommended configuration was arrived at. We devote much of the 
remainder of this chapter to the description of these tests and their results. For greater clarity, we 
depart from chronological order and describe the tests involving the time-domain and frequency- 
domain versions of ANC modules separately in the following two subsections. 
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4.2 Testing of Time-Domain Algorithms 

Referring to Table 4.1, test ANC 1 conducted initial tests of the basic segments of code 
comprising the time-domain, ARMA or ARMarkov-based neural replicators. Following success 
of these code segments on simple test cases, a SISO ARMA replicator was tried in ANC 2. For 
further reference, Figure 4.2-1 gives simplified block diagram conventions for description of the 
time-domain neural replicators. The simple block at the lower left of the figure represents an FIR 
time-domain replicator (given in expanded ANC notation at the lower right). Note that the 
adaptive speed in this block is replaced by unity so that in the diagrams that follow we can 
indicate the adaptive speed explicitly by multiplying the error signal by the appropriate factor. 
ANC 2 implements an ARMA-based replicator consisting of two such FIR blocks, with inputs 
equal to the disturbance system input and error sensor (one of the accelerometers) output and 
outputs summed to form the replicator output. Using a shaker at the MHPE base, good 
agreement was obtained with FRF data. With this success, ANC 3 extended the ANC 2 scheme 
to a 2 input/2 output MIMO replicator. After chasing down an unaccounted for time delay in the 
system, ANC succeeded in converging. However approximately 3 to 5 minutes was required for 
convergence. The reasons for the sluggish behavior were finally determined only after a 
diagnosis (as discussed below) of ANC 10 results. 

In any case, after ANC 3 it was judged time to add feedforward control to the identifier function 
in ANC 4. ANC 4 is shown in Figure 4.2-2. The upper portion of the figure shows an ARMA- 
type replicator where the replicator forward-path output has the form: 

y(k) = [^(^ ~m) + G m d(k - m) + P m u(k - ro)] 

y = error sensor output (4 2-1) 

d = disturbance input 
u = control actuator command 

The "algebraically constrained" ARMA controller is a deadbeat controller formed from the 
identified ARMA coefficient vectors F, G and P. Note, from the bottom half of Figure 4.2-2 that 
the control command, u , at time k is given by: 

u(k) = -y Fjr [F m y(k + ]-m) + G m d(k + l-nt)]+^ j P m u(k +1 - w)l (4.2-2) 
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Table 4.1: Candidate ANC Designs—Test Experience 



Configuration 

Result 

ANC 1 

Off-line checkout of delay line replicator 
code segments 


ANC 2 

SISO series-parallel ID 

Used a shaker input, got good match of I/O 
behavior using feedback system output. 

ANC 3 

2x2 MIMO series-parallel ID 

Required 3-5 minutes to converge. 

ANC 4 

2 x 1 MISO ID with algebraically 
constrained ARMA control 

Unstable - puzzling because of apparent good 
input/output behavior match prior to control. 
Note possibility of non-minimum phase 
model. Hard to debug (but see ANC 10 
below). 

ANC 5 

Off-line checkout of FFT/replicator code 
segments 


ANC 6 

SISO FFT/replicator ID 

Matched primary path well with a 1 Hz 
resolution. 

ANC 7 

2 x 1 MISO ID with FFT/replicator, 
control via division of identified primary 
and secondary path weights to form FFT 
based controller 

Distributed 5-30 Hz band, error and drift in 
the weights below 5 Hz caused large drifts in 
control => no appreciable attenuation overall. 

ANC 8 

2 x 1 MISO ID, control via pre-filtered 
replicator adaptation of FFT controller 
weights 

Drifting/unstable controller adaptation. 

ANC 9 

SISO FFT weight filter with fixed, user- 
selected FRF 

Verified ability to implement a specified FRF 
with this type of FFT filter. 

ANC 10 

2 x 1 MISO series-parallel ID, control via 
Explicit FRF computation and assignment 
to FFT weight controller 

Revealed serious roundoff concerns with 
weight updates and delay-line filter 
computations. 

ANC 11 

SISO ID of a single frequency bucket 
using cross and auto spectral density 
averaging based on recursive FFT 

Computes good results as compared with 
spectrum analyzer, but would like a little 
more speed. 

ANC 12 

SISO ID of multiple frequency buckets 
using ANC 11 algorithm, estimates 
accuracy of each bucket ID, adjusts FFT 
control of each bucket independently 

Good results, desire increased speed and 
decreased complexity. Also, improved rigor 
needed in ID (really need twice the resolution 
of control FFT in order to improve 
reliability). 

ANC 13 

Simple SISO FFT replicator with 
adjustable learning rate for ID (twice 
resolution of control), control algorithm 
requires no ID accuracy estimate 

Reliable convergence on center of channel, 
but poor performance in between. 

ANC 14 

Replicator with adjustable learning rate 
also used for control updates 

Reliable convergence, good overall 
performance compromise achieved over 
whole frequency band. 
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Symbols: 



or x(k) = [x(k),x(k-l),..,x(k-n+l)] T 



y(k) = h T (k) x (k) 
h(k+l) = h(k)+ x(k)z(k) 


Figure 4.2-1: Block Diagram Conventions for Time-Domain, ARMA-Based Neural 
Replicators 



Figure 4.2-2: ANC 4, Series-Parallel ID Using ARMA Model and Algebraically 
Constrained ARMA Control 


When this is substituted into (4.2-1), it is seen that y(k) = 0. Thus, once F, G and P converge to 
the actual ARMA coefficients of the system, one expects that y(k) converges to zero. As 
indicated by the switch shown at the lower left of Figure 4.2-2, a random signal is first injected to 
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the actuator command in order to perform system identification, then ujfis switched on to 
execute control. Although the identification step was apparently successful and good matching 
of input/output behavior was obtained, y(k) did not converge to zero. This puzzling behavior was 
thought to be the result of identification errors leading to a non-minimum phase model. 

However, as in ANC 3, the fundamental cause of the aberrant behavior was not properly 
diagnosed until the results of ANC 10 were understood. 

We attempted in ANC 10 one last try to implement a time-domain type replicator and to 
diagnose the results obtained in ANC 4. In the following we consider ANC 10 in detail because 
the results did lead to an understanding of the time-domain results and reveal a significant design 
issue with ARMA-based replicators. 


The algorithm tested in ANC 10 is shown in Figure 4.2-3. The replicator is a standard 
ARMarkov replicator with delay period L. The controller is synthesized from the ratio of the 
FFT (see the following section, equation (4.3-1)) of the forward path impulse response, G, to that 
of the secondary path, P. The FFT based filter for generating u had been extensively verified on 
earlier tests (described in the next section) so that any problems encountered could be attributed 
to the replicator. By focusing 


y(k) = y y(k) 



I Pi e-in 


Figure 4.2-3: ANC 10 Series/Parallel Plant I.D. Control via Explicit Computation of 
Cancellation FRF 

on the replicator operation and acquiring a wide variety of diagnostic outputs during the test, we 
were able to determine that both the weight updates and forward path signal computations 
suffered from serious roundoff errors that prevent controller convergence. Use of larger delay 
periods, L, in the ARMarkov model formulation quantitatively reduced the degree of error but 
did not qualitatively eliminate the phenomenon. The sensitivity of time-domain replicators (both 
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ARMA and ARMarkov) to roundoff error due to finite machine precision was then verified for 
the conditions of ANC 3 and ANC 4 and the same qualitative phenomena were duplicated on 
simple test circuits and correlated with analysis. 

The reasons for ARMA or ARMarkov replicator sensitivity to machine precision can be 
explained as follows. The problem shows up even in the forward path calculations because of 
inherent ill-conditioning of ARMA-type filters for high sample rates and systems composed of 
very lightly damped modes with frequencies well below the sample rate. To explain this, first 
note that a multi-mode structural model with forcing input xM and sample data sensor output 
y(k) has the ARMA-like representation: 

y(*) = £ yM 

v.(* + .) = [2^,6..-^] 

where ct m ,y m are modal influence coefficients c.() ands.( ) denote cos(.) and sin(.) y m 


represents the contribution of mode m to y and: 


0 m 

■ 2n t 



= m' h mode frequency 



= sample frequency 

(4.2-4) 

(P m 

■ exp [- Ti "ir 2,c ] 



= m' h mode damping ratio 



The above expressions can be manipulated to produce an overall ARMA model. To give an 
explicit example, consider two modes with a, = a, = a, y, = y, = 0 and (P,=(P,= 9 . We get the 
ARMA model: 

y(k +1) = [4? 2 (c.0, + c.0 2 ) 2 - 29 2 (l + 2c.6,c.e 2 )] y(k -1) 

+ [- 4? 3 (c.e, + c.e 2 ) (1 +2c.e,c.e 2 )+ 2 ^. 0 , + c.e 2 )] y(Jt - 2) 

+ [45 4 (c.e, + c.e 2 f-IF 4 ] y(k - 3) -2f 5 (c.0, + c.0 2 ) y(k - 4) 

+^(s.e,+s.e 2 )jc(jt-i) 

+ [2^ 2 (c.e, +c.e 2 )(5.e, +5.e 2 )-2^ 2 ^(e, +e : )]x(it-2) (4 . 2 - 5 ) 

+ [- 4P ? (c.0, + c.e 2 ) 5 .( 0 , +e 2 )+ip\s.e l + s.e 2 )] x(k - 3) 

+ 2'P 4 (c.0, + c.e 2 ) ( 5 . 0 , + 5.e 2 ) x(k - 4) 


36 


ANC Final Report 

Now consider what is actually computed at each time step when implementing the above as a 
digital filter. In particular we concentrate on the disturbance input terms (in*). When the 
disturbance frequency band is well below the sample rate and £2, is large, x(k-l-m) can be 
approximated by x(k-l) so that input terms are basically the sum of the ARMA coefficients of *. 
All of the significant coefficients are of the same order of magnitude (same exponent) but their 
sum is very small. Indeed, using small angle approximations for the sines and cosines, the net 
contribution of the disturbance terms to y(k+l) is approximately: 

5 ( e i + 0 2 )( l -^) 2 *(*-!) 

(G. +Q 2 ) ( £2, V „ „ (4.2-6) 

= io,t nH r-^: J 4 ’ 

Thus, if £2, is very large compared to £2, and £2 2 and T| is small (of the order 1 % or less) then the 
input terms in the ARMA model will sum to a very small number less than the last significant 
digit retained and the sum will be rounded off to zero. This predicts that the ARMA filter may 
work at a sufficiently low sample rate but produce zero output for nonzero input for larger 
sample rates. The same kind of phenomenon should occur as modal damping ratios are reduced. 
These predictions were entirely verified experimentally by simulating ARMA filters programmed 
into the DSP card by the outputs of analog filters set up to emulate three modes at 10, 20 and 30 
Hz with 1% damping. The ARMA filter implemented at 1000 Hz produced no discernible 
output. The modal damping ratios were also varied to produce the behavior predicted. 

In the general case, for ARMA models of the form: 

y(k) = X «,*(*-*)-S b,y{k-i) (4.2-7) 

i = l ,=1 

the transfer function: 

F(z) = (a 0 z n +... + a n ) / (z” +b,z n ~' +...+b n ) (4.2-8) 

has its zeros and poles clustered near z-1 for very light damping and high sample rate. As in the 
two mode case, the distance of the zeros from z=7 is proportional to the damping and inversely 
proportional to sample rate. This can generally be expected to produce significant roundoff 
errors in accumulating the input terms in the case of multi-mode lightly damped systems. In 
addition to the forward path signal computations, similar roundoff error problems were identified 
in the weight update computations. 

In summary, the limitations imposed by machine precision on successful implementation of time- 
domain, ARMA-based ANC algorithms was finally revealed by our combined theoretical and 
experimental investigation. In the light of these results we considered three options for further 
development. One option is to retain the time-domain ANC formulation but reduce the digital 
system sample rate used on the MHPE experiments. This could potentially circumvent the 
machine precision issue but would also significantly reduce the controller bandwidth. A second 
option is to again retain the time-domain formulation but somehow increase the processor word 
length. However, this would entail the procurement of much more expensive hardware that is 
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inconsistent with later transitions of the technology into commercial applications. The third 
option and the one selected is to implement and test the frequency-domain ANC algorithms. As 
discussed in the next section, tests of the basic frequency-domain modules, particularly in ANC 9 
showed no appreciable sensitivity to machine precision. Furthermore the frequency-domain 
ANC formulation meshes well with both classical and modem frequency-domain design insights. 


38 


ANC Final Report 


4.3 Test Results for Frequency’Domain ANC Algorithms 


Referring to Table 4.1, we initiated in ANC 5 experimental investigations of frequency-domain 
versions of the basic ANC modules. This sequence encompassed ANC 5-9 and ANC 11-13. 


Frequency-domain versions of the ANC components are formed from time-domain replicators by 
insertion of the finite-Fourier transform (FFT) operator at the inputs or the outputs. lix(k) is a 
discrete time sequence, define the FFT as: 


*[*(*)] A “I e-""* 1 * 

“ M mm 1 


x(k-m + l) 


(4.3-1) 


M = 2 n for some integer N 
r = 0,1,...»M-1 


The inverse relation is: 

x(k-m + l) = ^ e 3 r [x(k)] 

r-0 


(4.3-2.a) 


and, in particular: 

= [*(*)] 

r-0 


(4.3-2.b) 


The vector (; 0 [x(*)].7 1 [x(*], [*(*]) is obviously related to (x(k),x(k - l). x(k - M + ])) T 

by a unitary matrix operator. This operator can be inserted into the inputs to a time-domain 
replicator and the conditions for convergence remain essentially unaltered. However, the nature of 
the approximation changes. For example, if a system with output y and input x is FIR, then a time- 
domain replicator would be based upon the (forward path) dynamic model: 

y(*) = 2 K, x(k-m); M = 2 N (4.3-3.a) 


where the h m are the Markov parameters. But we can also use(4.3-2.a) to write this in the form: 

**>-£ »,»,[*(«] ,4 - 3 - 3b) 

r=0 

This translates into a replicator in which x is first submitted to an FFT operation, then input to a 
ganglion connected to the output neuron with weight vector [/7 0 ,..., H^]. Assuming the ganglia 
and synaptic connections are constructed according to established ANC rules, convergence is readily 
assured and (4.3-3.a) and (4.3-3.b) would seem to be entirely equivalent. However, in the frequency- 
domain version, the weights H r are now the values of the system frequency response function at the 

FFT (discrete time) frequency points 0 r = r—; r = 0,1,..., M -1. To show this, denote the - 
M 
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transform of x evaluated on the unit circle by Z e [x], where the transform and inverse transform 
relations are: 

(4.3-4) 

Z e [x] = X e e [o, 2 7t] 

k=~oo 


Then, the z-transform of the FFT of x is: 


z.[?,[4*)]]=H'„(e-e r )z,W, 9, A r 


2n 


(4.3-5) 




M(l-e^) 


is unity for £, = 0,2rt and vanishes at £ = u—;n = l,2,... A/-1. Hence, at 0 = r—, 

M A/ 

(4.3-5) gives Z e [.? r [x(jt)]] = Z e [x]. Now for a linear time invariant system: 

Z,[y]=W(9)Z,U] (4.3-6.a) 

But, taking the z transform of (4.3-3.b) and using (4.3-5) gives: 


r=0 


(4.3-6.b) 


Comparing the above expressions and using the properties of W M (o -0 r ) we see that: 

r = 0,1, •••, M - 1 


(4.3-7) 


Hence, a replicator based upon (4.3-3.b) is, in essence, a frequency analyzer. (4.3-3.b) represents the 
dynamic system by an interpolation function approximation to its frequency response. 

Referring again to Table 4-1, test ANC 5 was a checkout of the code for the FFT operation and 
the code needed to implement a basic replicator using the model (4.3-3.b). Following the code 
implementation and checkout, ANC 6 tested a SISO FFT replicator using the disturbance channel 
and one accelerometer on the MHPE tower. Clearly, the parameter M in (4.3-1) determines the 
number of frequency points at which the frequency response is matched. The expression for 0 f 
shows these points to be separated by ^ and this is the "frequency resolution" of the FFT-based 
replicator. In ANC 6, the frequency response was accurately matched by choosingA/ such that 
, = sample rate in Hertz) is 1 Hz. 
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With good progress on ANC 6, we moved on to test a MISO (Multi-input, single output) FFT 
replicator and a simple form of controller in ANC 7. To describe the ANC 7 algorithms. Figure 
4.3-1 first gives simplified block diagram conventions. The FFT block illustrated at the top of 
the figure computes the FFT of the input x(k) using the fast, simple recursion relation: 


5,[x(k + l)] = e*' ?,[*(*)] +— [*(* + l)-;c(Jt-M + l)] 
M 

n a 2U 

0 r A r — 
r = M 


which is readily obtained from (4.3-1). The output, x(k), contains the real and imaginary parts, 
given alternately, for all the FFT components. 

The "divided" weights" device, shown in the lower half of Figure 4.3-1 synthesizes an FIR filter 
based on representation (4.3-3.b) where the frequency response values at the FFT frequency 
points (the Hr's ) are computed from the quotient of two given frequency responses. 

Obviously, both devices in Figure 4.3-1 are specialized types of linear filters and involve no 
adaptive adjustment of coefficients. 

With the above conventions, Figure 4.3-2 shows the algorithms tested in ANC 7. The top half of 
the figure shows an FFT version of an FIR replicator. In the frequency-domain, this model is 
based on the system model: 

Z,[y]= G(9 )Z, [d] + P(6 )Z, [«] (4.3-9) 


y = sensor output 

d = disturbance signal measurement 
u = control command to actuator 

G(0) and P(Q) are termed, in the acoustics literature, the "primary path" and "secondary path" 
dynamics and we use the same terminology here. The replicator in Figure 4.3-2 attempts to 
approximate Z e [y] at the M FFT frequency points. The z-transform of y is given by: 

Z »[>]=S G r w u (o-e,)z,ld] (4.3-10) 

+X W»(e-e r )z.[«] 
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Recursive FFT: - 


FFT(n,Pi,p 2 ) 


-► => *«(k) All x(k-r) e j £ U 

n r=0 

Re{xp|(k))] 0<pj<p 2 < ii 
Imgpjfk)} 


Re(iP2*)l 

ImGp 2 (k))| 


Divided Weights: 



y(k)=H(k) T x(K), 

(H2i>l) + j(H2i ) = (G 2i-l) + j(G 2j) 
(P2i-l)+j(P2i) 


i=l:Length(G) 


Figure 4.3-1: Some Block Diagram Conventions Relating to Frequency-Domain Neural 
Replicators 



Results/Conclusions: 

• fsample = 500 Hz, n=500, pj=0, P2=40, nswitch/f sam pi e = 180 sec 

• disturbance from 5-30 Hz 

• large low frequency drifts observed when control switched on, which used up dynamic 
range 

Figure 4.3-2: ANC 7, FFT Replicator for Plant ID, Control via Division of Identified 
Transfer Function Weights 
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Thus, if G r and P r converge to G(0 r ) and P(0 r ), respectively, then Z e [v] exactly matches the 
system frequency response at the FFT frequency points. Once replicator convergence is obtained 
and one has a good frequency-domain model, the optimal feedforward controller is such that: 


Z,[u]= 


G(e) 
p(. 0 ) 


z,[d] 


(4.3-11) 


The implementation of this is shown in the lower half of Figure 4.3-2, using the "divided 
weights" filter defined in Figure 4.3-1. With the values of the parameters indicated at the bottom 
of Figure 4.3-2, persistent errors and drift were observed in the low frequency weights (<5 Hz) of 
the replicator and, in consequence, large drifts were obtained in the control input. 

At first, it was thought that the above might be due to the use of the "divided weights" device for 
computing the control and due to the representation of the HR by an FIR replicator. Hence the 
algorithms of ANC 8, shown in Figure 4.3-3, were tried. In this case, the replicator is an ARMA 
model and the controller is an FFT version of an HR replicator, which computes the needed 
cancelling control input given the identified secondary path dynamics,/’. 

Basically the same drift in the low frequency weights was observed in ANC 8 as in ANC 7. It 
was concluded that the use of the recursive control computation, shown at the bottom of Figure 
4.3-3, offered little benefit and that the use of an FIR versus HR replicator was irrelevant to the 
difficulties experienced in ANC 7. In any case, a significant portion of the low frequency drift 
was traced to the round-off problems with ARMA replicators discovered in ANC 10. 

Our next hypothesis was that the drift problems experienced in ANC 7 were due to the fact that 
replicators based on (4.3-3.b) are only " half frequency domain". Replicator inputs are FFT 
components but their outputs are still representations of the time domain output of the system 
being identified. Implicit in this setup is that all FFT components and corresponding frequency- 
domain weights have to be employed—including the lowest frequency components. As the low 
frequency components are most subject to measurement errors and drift and are associated with 
the slowest convergence rates anyway, it is, in retrospect, no surprise that the replicators could 
not approximate the low frequency behavior within reasonable learning periods. However, to 
achieve good feedforward control over the frequency band of the disturbances, it is not actually 
necessary to approximate the system transfer functions everywhere in frequency but only in the 
disturbance band, starting at some reasonable low frequency cutoff. Thus, our response to the 
difficulties of ANC 7 was to devise replicators that only operate on a designated frequency 
band—thus eliminating the processing of very low frequency components. To do this, it is 
necessary to use fully frequency domain replicators for which the inputs are FFT componentsa/jd 
the output is a representation of the system frequency response (not its time-domain response). 
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Results: 

• unstable learning, but I.D’d plant seemed good. (!?) 

• f s = 500 Hz, n=15, N=500, d=0.01 °= H = 0.01, d(k) 5-30 Hz random 120 seconds of 
learning 

Figure 4.3-3: ANC 8 2x1 Series Parallel I.D. XLMS FFT Weight Adjustment 

In place of the model (4.3-3.b) suppose we approximate [>’(&)] - H 6r J,. [*(&)]—i.e. attempt to 
replicate an FFT component of y(k) rather than y(k) itself. The ANC replicator corresponding to 
this setup is shown in Figure 4.3-4.a, where < • > denotes the time average: 

It,' . (4.3-12) 

< *(*)> = — 2, x(k-m) 

1 m=o 


and we assume that the actual system response is given by (4.3-6.a) 

Following the rules of ANC construction, the update rule for H %r in Figure 4.3-4a is given by: 

H„(n + 1) = H„in) + ° U[y]K[ U ]- H„ |^[>]| ! ) (4 ' 3 ' 1 

(KMI / 

and for a « 1 we may use the averaging approximation: 


44 


















ANC Final Report 


" H„(n +1)» H„M + K Kl'f ] 

f(?-[y]r[4 


(4.3-14) 


= H 6r (n) + a\ 


W 




Now consider the quotient within brackets. By Parseval’s Theorem and (4.3-5): 

m [l\y\^[u\) = 1 ^ 

T - 


^iim (?,[y]CM)=^f de|^(e-e^//(e)|z 6 M 2 | 


(4.3-15) 


Since W M {^} assumes nonnegligible values only near £ = 0 and supposing adequate frequency 

resolution (— sufficiently small), we obtain: 

M 


Wy]VM> =/ V>(W"f) 


(4.3-16) 


Thus the update equation becomes: 

H„<» +1) = H.M + o(H(9,)- <4 3 ' I7) 

and with a < 2, H 6r converges to //(0 r ). Thus the replicator of Figure 4.3-4a directly estimates 
the frequency response at one FFT frequency. Notice that the replicator operates on only one 
FFT component. As many or as few frequency components as are desired can be processed in 
parallel by a multiplicity of such mutually independent replicators. Besides rapid convergence 
and simplicity, this scheme allows one to identify systems over a specified frequency window 
and exclude unnecessary low frequency components. 

Figure 4.3-4b shows the forward path signal flow for combining several "frequency bucket" 
replicators, and effecting the inverse FFT to synthesize a time-domain output. This is needed 
later to synthesize the proper control inputs. As a preliminary step, the forward path algorithm of 
the "FFT weight filter" in Figure 4.3-4b was implemented and successfully tested in ANC 9. 

Subsequently, ANC 11 tested the capability to identify the disturbance to error sensor transfer 
function value at an FFT frequency point (in a single "frequency bucket") preparatory to a 
complete identification exercise in ANC 12. 

For greater stability and accuracy in ANC 11 we actually employed a "batch" version of the 
algorithm in Figure 4.3-4a. Figure 4.3-5 shows the system tested in ANC 11. The algorithm is 
depicted in previously established ANC block-diagram notation and is devoted to identifying the 
transfer function at a single frequency point. Identification over a frequency band is 
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^ r [x]- 


-» q> 4- 


-■^y] 


-4-CW- 


(i = a/< |^ r [u ]| 2 >, a = j«l 

(a) 


/ M 


: . 1 

diag {HG^} 

\ y N w " 1 

k 


> m -►y(k) 


e ^ (1,1,...,1) T 
(b) 


Figure 4.3-4: Devices for Fully Frequency-Domain ANC Components, (a) Single Channel 
Replicator and (b) Forward Path Dynamics Combining all Channels to Produce a Time- 
Domain Output 



Note: AH 6r =\ix H e, p = cc/|jx|| 2 , a<2 

Exit if: lie’ll 2 / (( m - l)||x|| 2 ) < X 2 (= required mean-square accuracy of H 6r ) 

Figure 4.3-5: Batch Version of Fully Frequency-Domain ANC Replicator Tested in 
ANC 11 
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accomplished by an array of modules of the type shown, all running in parallel. The basic inputs 
to the replicator assigned to identify at frequency d r are the FFT components of the input and 
output; J r [.x(k)] and £,[>(*)], respectively. The input and output ganglia are connected by a 
synaptic bundle with weight matrix constrained to be proportional to the identity matrix. 

Because of this constraint, the weight increment is proportional to the inner product x H "e as 
shown in the formula in Figure 4.3-5. This frequency-domain identifier basically time averages 
the results that would be given by the algorithm of Figure 4.3-4a, thereby achieving good noise 
filtration properties. The number of time samples used in the averaging is given by the parameter 
m. Assuming the learning rate constraint, a, is less than 2, H 6 will converge to H (0 = 0 r ) 
supposing that x and y are related as in (4.3-6.a)). In the subsequent work, a=l .0 was used 
throughout. Because the weight matrix is diagonal, the computational burden is actually rather 
small and, as will be seen below, a large number of frequency channels of the form of Figure 4.3- 
5 can be implemented simultaneously on the DSP card. Finally, as indicated at the bottom of 
Figure 4.3-5, we use x and e to compute an estimate of the mean square identification error. 

This is compared with X 2 , the user-supplied requirement for mean square accuracy. If this test is 
passed, the algorithm terminates and announces the determination of a frequency response 
estimate of the desired accuracy. 

The above algorithm was tested in ANC 11 and the accuracy of the results was found to be very 
good when compared with the results of the spectrum analyzer. The algorithm also showed its 
ability to monitor its own accuracy. With this encouragement, ANC 12 went on to demonstrate 
complete identification over a frequency band. In ANC 12, the algorithm of Figure 4.3-5 is 
simply repeated for every FFT frequency point within the desired frequency band, and the results 
for H m ,..., H 6n are assembled to form the FIR filter representing the identified system. 

As noted in Table 4-1, results for ANC 12 were quite good. For example, using 500 delay stages 
in the FFT operation with a sample rate of 250 Hz (giving a 0.5 Hz frequency resolution), the 
transfer function from the disturbance command to one of the MHPE tower accelerometers was 
identified over the 5 to 25 Hz band with 51 "frequency buckets." 

Figure 4.3-6a, b compares the transfer function gain and phase plots of the identified filter 
(dashed line) with FRF data (solid lines). The match is excellent for both gain and phase. In 
neither ANC 12 nor in any subsequent tests were roundoff error problems observed. 

With successful performance of the frequency-domain system identifiers, we proceeded to 
implement the frequency-domain control algorithms. Figure 4.3-7 shows the overall controller 
structure in which each frequency channel is identified and controlled separately and in parallel 
with all other channels. First, y and d are subjected to FFT operations. The FFT blocks in Figure 
4.3-7 show the digital filter realization of the recurrence relation (4.3-8). The 
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ANC12 Results: 500 taps at 250Hz => 51 buckets 



Frequency (Hz) 


Figure 4.3-6a: Gain of Disturbance to Sensor Transfer Function. Frequency-Domain ANC 
Replicator (dashed line) is Compared with FRF Data (solid line) 












Phase (Degrees) 
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Figure 4.3-6b: Same Comparison as in Figure 4-6a Except that the Transfer Function 
Phase is Shown for ANC Replicator (dashed line) and FRF Data (solid line) 

r th FFT components of y and d are denoted by y r and d r , respectively. Each pair of FFT 
components having the same index (the same center frequency) are handled separately. For 
r = 1,..., M , each pair (y r d r ) are input to the rth channel subsystem consisting of an rth 
channel identifier in series with an rth channel control adaptor. The channel identifier outputs an 
estimate, G r of G r (0 = dr) where G is the closed-loop transfer function matrix from d to y and 
also passes along d r . With the inputs d r and G r , the channel r control adaptor outputs its 
component, u r of the actuator drive signal vector. Finally, outputs of all the channel control 
adaptors are summed to form the total actuator command signal vector, u(n). 
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Figure 4.3-7: Overall Controller Processes the Separate Frequency Channels in Parallel 
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To understand the details of the identifier and control adaptor devices in Figure 4.3-7 we sketch 
the control issues that must be addressed. First, repeating (4.3-9), the system being controlled 
has the frequency-domain form: 

Z,[y]= G(9)Z,[d]+ />(0 )Z, [«] (4.3-18) 

y€lR/, weIR m 

The feedforward controller is a linear filter with input d of the form of (4.3-3.b): 

«(»)=!>,?,[</(»)] <4 ' 3 ' l9) 

k=0 


where H k 3 t [d(n)] corresponds to u k in Figure (4.3-7). Substitution of (4.3-19) into (4.3-18) 
and use of (4.3-5) gives the closed-loop relation: 

Z„[y]=G(e)Z,[d] 


G(e)= G(e)+ p(9) X h , M'„(e -e,) 

*= o 


(4.3-20) 


so that G defines the closed-loop transfer function matrix. From the above, it is seen that the 
basic control strategy is to (somehow) identify the quantity P~'(d r )G(B r ) for all dr in the desired 
control band, then set H r =- P~‘(0 r )G(6 r ). This makes G(0) pass through zero at each FFT 
frequency point. With adequate resolution, this strategy produces the desired level of response 
suppression over the required frequency band. 

We desire to carry out the above control strategy so that the identifications are performed and the 
cancelling control is synthesized simultaneously without the use of extraneous test signals (e.g. 
broadband dither introduced into u to identify P(9)) and without interruption of normal 
operation. There are two basic approaches: 

(A) Using only d k , u k and y k , identify g( 0 t ) and p(d k ) separately, then generate 

H k =-p-‘(e i )G(eJ 

or: 


(B) Identify 2 k A P"'(0 A .) and g( 0*) and then generate H k from: 

H k (n + l) = H k (n)~ £. k G t (n) (4.3-21) 

The recursive scheme (4.3-21), in effect, identifies - P~'(p k )G(d k ) to produce the desired 
controller. To see this, use the expression for G in (4.3-20): 
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. H k (n +1) = H k (it) -S t [g(0,)+ Pfa )H k («)] 

= H k ( n) - [p~ ] (e 4 ) G(e k )+ H k («)] (because & k = P~' ( 0 ,)) 
=-p(eJc(eJ 

Approach A has had considerable success, especially in tonal cancellation. However this scheme 
has a "slow mode" of convergence that, over extended operation, progressively slows down the 
convergence rate. On the other hand, approach B does not exhibit this behavior and is very 
robust because (4.3-21) can converge to the cancelling controller even when the error 
- P~ ] (0, | is quite large. Thus, approach B was chosen for the present study. 

The basic, preliminary identification task needed to implement approach B is the identification of 
G (0,), r = 1,..., M. Each of the G (0 r ) is determined by one of the channel identifiers in 
Figure 4.3-7. The channel r identifier has the form shown in Figure 4.3-5 with d r and y r as 
inputs and G, as the output estimate of G (0 r ). For this "batch" frequency-domain identifier, it is 
most convenient to run the algorithm in separate time blocks, each of length A/. In other words, 
the estimate G t is held constant for M time steps after which a new estimate is produced. If M is 
roughly the correlation time of each FFT coefficient, the algorithm of Figure 4.3-5 takes L steps 
to arrive at a converged estimate where L is typically 3 to 5. With this time-block mode of 
operation, the control adaptors update the H r 's only once every "identification epoch" of LM 
time steps. Thus, in describing the control adaptor operation, the time index V is used to 
characterize the values of 2 k ,H k and G k on the n identification epoch. 

Assuming the channel identifier has arrived at an acceptable estimate of g( 9,) at time "n", 
we have from (4.3-20): 

U k (n) = G(0*)+ P(Q k )H k (n) +T\ k (#i) (4.3-22) 

where T|, represents sensor errors, instrumentation noise and identification error. A bound on the 
noise fluctuation can be assumed in the form: 

max |n(«)-Tl(« -1)|| < cr (4.3-23) 

n 

Since is largely the consequence of the electronics hardware being used, we may estimate (on 
the basis of test results) a value of c. This was the procedure used in early tests; subsequently the 
algorithms were augmented by an on-line estimation of a. Also, for future reference, define the 
notation: 

A»(n) A •(n)-»(/i-l) (4.3-24) 

and note from (4.3-22) that: 

A U k (n) = P(e k ) A H k + At), (ii) ( 4 - 3 ' 25 > 

With the above relations, we can describe the operation of the k th channel control adaptor, which 
is shown in Figure 4.3-8. We use the established ANC block diagram conventions. The main 


52 




ANC Final Report 


input is G k (n). A ML means a time delay back to the previous identification period. The lower 
portion of Figure 4.3-8 is a device that makes use of the correlation between fluctuations of G k 
and fluctuations of H k to determine P~' (0 t ). Essentially, this is an identifier for a linear 
mapping from A H k to A G t (n). From (4.3-25) it can be seen that 3 t converges to P~' (e k ) 
(since, ignoring noise, A H k = P~ ] (0* )A G k ). Given G k and a copy of 2. k , the upper half of 
Figure 4.3-8 implements (4.3-21). Specifically, following the block diagram conventions; the H k 
update is: 

H k (n + 1) = H k (n) - 3 k (n + \)t k (n) 

—i.e. H k is updated just after 2 k is. Note that the learning rate for the 3 k adaptation is adjusted 
during operation according to the formulae as given at the bottom of Figure 4.3-8. This prevents 
undue noise sensitivity when convergence is nearly complete and G is small. Here, o is the 
bound on noise fluctuation, (4.3-23) and € is a prespecified tolerance on the P~' (0*) 
identification, i.e.: 

€ = required upper bound on ||l - P(0 k )2* (n)|| (4.3-26) 

For the above algorithm, we can state the following result: 


Suppose d(k) is a stationary random time series and P(z ) is stable, 
strictly proper and minimum phase. Then for the algorithm of Figure 4.3-8 and 
for any initial values of $ k and H k : 

lim ||/ - P(0 t ) 2 t (n)|| <e 
n —> 

lim ||(7 t («)| < a / (1- €) 


a. (4.3-27) 

b. (4.3-27) 


with probability one. 


Thus the algorithm is guaranteed to achieve the desired accuracy in identifying the "inverse 
secondary path", P~ ] (6 k ) and in reducing the closed-loop transfer function down to the 
instrumentation noise floor. 


The identification and control algorithms described above were implemented in ANC 12. This 
test sequence concentrated on a single control loop involving one of the MHPE tower-mounted 
LPACTs and its colocated accelerometer. Some six tests were performed and the results were 
generally quite good. For example, Figure 4.3-9 shows one case in which the desired control 
range is 5 to 20 Hz and the frequency resolution (interval between FFT points) is 0.5 Hz with a 
sample rate of 250 Hz. It is seen that the dominant group of modes near 10 Hz is reduced by 
approximately 20 dB. The result shown in Figure 4.3-9 was obtained after 3 iterations of the 
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G k eC" 

u k e C m , m > (- n 
d k e 



In the H k Adaptation, the Learning Rate is: 

1, |a (j, (n)| > a/ e, e< 1 
oc = * * * 

0, otherwise 

a = bound on noise fluctuation 
e = ID error tolerance 

Figure 4.3-8: Selected Structure for the k^ 1 Channel Control Adaptor is Essentially an 
Identifier for P’ ] (e, ) and P~ l (e, )g( 0 t ) 

algorithm in Figure 4.3-8. In this scheme, most of the computation time is taken up by the 
identification process since the length of each identification epoch must be more than the 


54 















ANC Final Report 


correlation time of each FFT component ofy or d. Since the correlation time is roughly the 
reciprocal of the frequency resolution the overall convergence time is inversely proportional to 
frequency resolution. 

Despite favorable results on ANC 12, the closed-loop response after convergence was found to 
vary noticeably from run to run. After further analysis this was traced to two causes. First on 
each epoch, the identification algorithm in Figure 4.3-5 was run for a prespecified time interval 
and a fixed number of epochs was used for each control update. Thus the identifier provided an 
estimate of G k to the controller without regard to the actual accuracy of the result. Secondly, in 
the above scheme the controller adjusts its gain for the same number of FFT components as are 
included in the estimator—i.e. the frequency resolutions of the identifier and controller are the 
same. In this case, there is discernible interaction among the frequency channels, particularly 
where G(0) varies rapidly over a frequency bucket. Such interactions cause nonrepeatable 
variations in the final controller. To reduce these interactions we observe that if the plant is 
adequately approximated by an FIR filter of M delays and the controller is also an FIR filter of 
the same order then the identifier model of the close-loop must have 2 M delays—i.e. the 
frequency resolution of the identifier ought to be twice that of the controller. When this 
condition imposed, it is clear that cross-channel interactions are minimized. 

Therefore in ANC 13, we modified the algorithm of ANC 12 by (1) running the identifier so that 
it terminates only when the mean-square error estimate falls below X 2 and (2) setting up the 
controller so that it addresses every other FFT point that is processed by the identifier. 
Modification (1) ensures estimates of G(0) of uniform quality. Modification (2) gives the 
identifier twice the frequency resolution as the controller. With these improvements, ANC 13 
converged somewhat more quickly and far more reliability than did ANC 12. The algorithm 
development would have concluded here if it were not for one fundamental difficulty. 

The problem was that because each frequency channel was treated as if it were independent and 
the fully frequency-domain ANC replicator (Figure 4.3-5) was used, the ANC 13 algorithm 
tended to control the system at the centers of the frequency channel but not elsewhere. What 
often occurred was that deep (>20 dB) nulls were created at the center of the various channels, 
but elsewhere the closed-loop response was near the open-loop levels. 

One of the essential difficulties was traced to the fact that the fully frequency domain replicator, 
Figures 4.3-4;5, when used to identify the closed-loop transfer function G(0) in equation (4.3- 
20), actually introduces coupling among the several frequency channels. Recall that the 
replicator works by correlating £ r [y] with £ r [</]. Using (4.3-5) and (4.3-20), we find: 

r r n / x T M-l / xl (4.3-28) 

z e W M] = w 2 m ( 0 - Q r) |G(0) + m 1 q H k W M (e -ejj z 0 [d] 


From this, it is evident that -5 r [v] has significant frequency content outside the t*h frequency 
channel coming from the k^ (k*r) control channel input. This is still a source of significant 
variability in the identification and control results. 
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A secondhand more significant difficulty is that the frequency domain replicator and independent 
channel controllers take account of the closed-loop response only at the channel center 
frequencies. The algorithm has no mechanism to recognize closed-loop response in between the 
center frequencies. This primarily accounts for the qualitative results noted above for ANC 13. 

Upon reconsideration of the various replicator forms tested, it was decided to use the partial 
frequency domain versions of ANC of the type shown in Figure 4.3-2. These are based on the 
model form (4.3-3.b) and employ FFT inputs but time domain outputs. Figure 4.3-10 shows the 
mixed time-and frequency-domain algorithm used for ANC 14. This constitutes the finally 
successful algorithm adapted as the basis for our final MHPE demonstrations. 

The algorithm shown in Figure 4.3-10 simultaneously identifies and controls the system. As in 
ANC 12 and 13, closed-loop system identification and control adaptation occur in epochs which 
are set to a fixed period (4 minutes in most experiments). G(0) is identified (see the bottom of 
Figure 4.3-10) using a model of the form (4.3-3.b). The input is the 2n FFT of the disturbance 
measurement while the output is an attempted replica of the error measurement y. As discussed 
in connection with (4.3-3.b), when G converges, it converges to the closed-loop transfer function 
from d to y at the FFT frequency points. Note that the time-domain control signal, u(k ), is 
synthesized in a similar way—i.e. using an n point FFT of d. Thus, as in ANC 13 the frequency 
resolution of identification is twice that of control. 

Using the difference of G over different identification epochs and the current value of the control 
gain vector, H eC", the secondary path dynamics transfer function ,P is synthesized. These 
computations rigorously produce the required secondary path response, assuming the closed-loop 
replicator converges perfectly. This approach is also relatively insensitive to noise due to the fact 
that the differences in G and the control transfer function (i.e. H) magnitudes are quite large 
compared to the noise, especially after several adaptations have occurred. 

After subjecting the output of the P T block to an n-point FFT operation, the result is input to 
another replicator with synaptic weight vector dH T . This replicator matches its output to the 
time series signal -y(k). In other words dH converges to the change in the control gain needed to 
make dH{Q)P(Q) equal to — G(0). Thus, when dH is added to H, the desired feedforward 
controller is obtained. This is done during an identification epoch while maintaining the on-line 
gain, H, fixed. When identification epoch £ is complete, dH f - is added to H - and H is updated. 
This approach to control update ensures that the time series y(k) converges to zero in a least 
square error sense. Hence, the algorithm attempts to minimize the closed-loop frequency 
response magnitude both in the vicinity of the FFT points and between these points. Thus, the 
fundamental difficulty of ANC 13 is avoided. By using the model form (4.3-3.b), the algorithm 
also avoids the sensitivity to machine precision observed for ARMA models. Finally, we note 
that the controller gains are updated at the beginning of each epoch with a perturbation which is 
computed off-line. This guarantees that the control and identification processes are completely 
independent and produce well-defined estimates of the desired quantities. 

The ANC 14 algorithm was the basis for the final experimental demonstration. We review the 
results achieved by this algorithm in the next section. 
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y(k) 


• Superscript £ (e.g. ) refers to learning epoch £ 

• Learning for dH and G run simultaneously, and are reset for each epoch 

• H is fixed during learning epochs 


Figure 4.3-10: Mixed Time-and Frequency-Domain Algorithm for ANC 14 

4.4 Experimental Results for ANC 14 Algorithm-Final MHPE Demonstration 

The final results for ANC continuous spectrum control are presented in this section. These 
experiments use the ANC 14 algorithm described above. The ANC 14 algorithm was run using 
several actuator and sensor configurations on the MHPE testbed. Results were found to be 
affected by linearity, signal-to-noise ratio, and the inherent limitations of the filters used to 
control the system. 

All experiments were run with a 250 Hz digital sampling frequency and learning epochs of 4 
minutes (240 seconds). In what follows, the parameters of interest are df (the control resolution), 
fj (the lowest control frequency), and f2 (the highest control frequency). Note that cancellation 
can only be expected after three epochs (12 min) of computation, since the first epoch is used to 
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compute open-loop response, the second is used to compute the first secondary path estimate 
with an arbitrary control input, and the third is used to compute the first control perturbation. 

In one experimental configuration, the disturbance was injected into the D/A channel of LPACT 
#1 located on the MHPE tower while sensing and control were executed using LPACT #2 on the 
tower. For good signal-to-noise ratio, a disturbance level of 500 mV was used. 

Results for the case of df=0.4, fj = 3.5 Hz and f2=24 Hz are shown in Figures 4.4-1 through 
4.4-3. Attenuation of 13.4 dB rms was obtained after 20 min of adaptation, with approximately 

23.5 dB obtained on the major structural mode at around 16 Hz. 

At this point, it was decided to try a direct comparison with the earlier ANC results. Thus, in the 
second test configuration investigated here, the disturbance was injected into an LPACT mounted 
on the main petal of the MHPE, as before on ANC 12. The disturbance LPACT had a much 
lower force output than the shaker, so an input level of 1 volt rms was used over the same 7.5-20 
Hz range as previously. This still did not give a very high signal-to-noise ratio for identification, 
however. In order to remedy this deficiency, it was decided to increase the gain of the LPACT 
A/D output by a factor of ten. The final set of experiments was then run to establish the best 
achievable performance. The parameters were df=0.4, fj=3.5 Hz and f2=24 Hz. Figures 4.4-4 
through 4.4-7 show the results, with perhaps the most impressive case being Figure 4.4-7, the 
result after 24 minutes of learning. Attenuation overall was 14.1 dB (probably more if noise 
spikes were averaged out), with approximately 23 dB attained on the primary resonance at 

14.6 Hz. 

In summary, the ANC 14 algorithm appears to attain reliable performance limited only by 
theoretical implementation considerations and the signal-to-noise ratio. 
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Petal to LP1, df_ctl = 0.4, 3.5-24 Hz CTL, lOx A/D, 16 min Learning 



Figure 4.4-5: Same Case as in Figure 4.4-4. Open-Versus Closed-Loop Response After 
16 Minutes of Adaptation 
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Petal to LP1, df_ctl = 0.4, 3.5-24 Hz CTL, lOx A/D, 20 min Learning 



frequency (Hz) 


Figure 4.4-6: Same Case as in Figure 4.4-4. Open-Versus Closed-Loop Response After 
20 Minutes of Adaptation 
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Petal to LP1, df_ctl = 0.4, 3.5-24 Hz CTL, lOx A/D, 24 min Learning 



Figure 4.4-7: Same Case as in Figure 4.4-4. Open-Versus Closed-Loop Response After 
24 Minutes of Adaptation 
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5. Amended Option 1 Program Overview 

The Amended Option 1 Program addition was dedicated to the development, delivery, and 
demonstration of an ANC system capable of canceling three tonal vibrations in a fault-tolerant 
fashion using the three ACESA struts on the ASTREX testbed at the Air Force Phillips Lab. The 
contract addition was made at the end of June, 1995 and extended the contract to the middle of 
September, 1996. 

This multi-tone demonstration system had two important features that had not been demonstrated 
previously on the ANC program: 

1. Multiple sensors (3) and actuators (6) were used. The sensors were fewer in number than the 
actuators and were not collocated. 

2. In addition to controlling the ASTREX structure without any prior modeling data, the system 
was designed to be capable of recovering from actuator failures autonomously. 

The system that was developed used an architecture that combined analog and digital ANC 
elements. This resulted in very low throughput required—the sample rate of the digital computer 
was set at 1/20 Hz. The hardware that was used to implement the algorithm was vastly more 
capable than required to execute the computations needed by the ANC system. It was chosen to 
minimize development cost because of commonality with another program. 

The system was demonstrated during two test visits, one in January and the other in June of 
1996, each for a week. 

In the following, we describe the demonstration system in Chapter 6. We then go into detail 
describing in mathematical terms the algorithms tested in Chapter 7. Chapter 8 describes the 
experiments that were run and the results achieved. Chapter 9 describes the software that 
comprised the demonstration system. Chapter 10 gives an overview of the hardware, focusing in 
some detail on the custom analog cards that were delivered. Finally, Chapter 11 gives a 
summary of the conclusions and future directions of research for the whole ANC program. 
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6. ASTREX Demonstration Description 

In this chapter we describe the ASTREX demonstration architecture in some detail. We begin 
with an overview and then proceed to the subsystems that make up the whole. 

6.1 Overview of the Demonstration 

A summary block diagram of the demonstration is shown in Figure 6-1 below. It consisted of six 
major components: the ASTREX structure, the ACESA secondary mounting struts and 
electronics, the Linear Precision Actuators (LPACTS) and servoamps, the three tone disturbance 
source, the motion sensors and amplifiers, and the ANC control system. 
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Three sinusoidal disturbances are generated, scaled, and summed in the disturbance source block. 
These sums are fed to two of the three LPACTs on the secondary tower of the ASTREX test bed 
as system disturbances (one of the three LPACTs was not functional at the time of test). 
Meanwhile, three pulse trains, synchronized with each of the three sinusoids, are passed to the 
ANC controller. Thus, the controller knows the frequency, but not the amplitude or phase of the 
disturbance inputs to the LPACTs. The controller synthesizes six canceling control signals to 
feed to the ACESA struts, one for each bending axis of the three struts. The resulting motion of 
the structure, which is to be canceled out, is sensed and fed to the ANC controller as a training 
signal. 

The demonstration emulates a precision optical system that is disturbed by independent 
disturbances such as might originate from three control moment gyros (CMGs) with slight 
imbalances. Motion of key optical subsystems that contribute to line-of-sight (LOS) error are 
sensed and canceled by the ANC controller using actuators capable of controlling the relevant 
substructures (in this case, the ACESA struts that support the secondary mirror assembly). 

In addition to the components shown in the Figure, a separate system was installed prior to the 
second test visit to enable visual observation of the cancellation system. This system is depicted 
in Figure 6-2 below. 



Figure 6-2. Performance Visualization System for the ASTREX Demonstration 


Two laser sources are mounted near the cylindrical telescope mass simulators in the figure. 
These aim at two small mirrors mounted on the secondary assembly. The laser beams are 
reflected on to paper targets mounted on the primary structure as shown in the figure. This 
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provides complete monitoring of the effects of the secondary mirror assembly on the LOS 
performance, and is entirely independent of the sensors used for training the ANC. 

We now discuss the various components of the demonstration in somewhat greater detail. 

6.2 ACESA Struts and Electronics 

The ACESA struts were developed by TRW as part of an earlier effort. The interface to the 
ANC system is via six analog inputs, which command the three struts, each in two bending axes. 
The struts are labeled for convenience the 12:00 strut, the 4:00 strut, and the 8:00 strut according 
to their positions as viewed by an observer looking in the -Z direction (down the boresight from 
the secondary mirror assembly towards the primary). 

Each strut has two inputs: one to cause motion of the secondary in the X direction, and one to 
cause motion in the Y direction. For safety reasons, the inputs were limited to ±3 V to prevent 
damage to the structure or struts. 

6.3 Motion Sensing for Adaptation 

Two different schemes were used during testing to measure the motion of the secondary mirror 
assembly. The first scheme consisted of two accelerometers mounted in the X and Y directions 
in the vicinity of LPACT 2. This scheme was sensitive to lateral motion of the secondary and 
also torsional rotation about the Z axis. It was used during the first test visit and part of the 
second. Two different sensitivities were used, 10 mV/g and 100 mV/g. 

The second scheme consisted of three angular velocity sensors mounted near the mirrors on the 
secondary. These were insensitive to translation of the tower and provided a much better 
indication of the effect of the secondary on the visually observed spot motion of the laser beams 
on the targets. The velocity sensors had sensitivities of approximately 400 mV-sec/rad. 
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6.4 LPACTS and Disturbance Sources 

The disturbance generation system works according to the block diagram shown in Figure 6-3 
below. 


synch synch synch 

signal 1 signal 2 signal 3 



Figure 6-3. The Disturbance System for the ASTREX Demonstration 


Three signal generators provide individually adjustable sinusoids, which are summed and input 
to the two LPACTs used for disturbance generation. The summing function is provided by a 
simple op-amp circuit. The same disturbance signal is input to both LPACT servo-amps. As 
mentioned previously, LPACT 1 was inoperative due to an intermittent connection. 

The LPACTs provide approximately 10 lb. maximum output for a maximum input of 
approximately 5 volts at 10 Hz. Slightly more output is available at higher frequencies. 
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6.5 ANC System 

The Adaptive Neural Control System has the structure shown in the block diagram of Figure 6-4 
below. Except for the VME computer peripherals (disk units and terminal), the whole assembly 
fits in a VME card cage. 

This implementation of the ANC system is based on the observation that if the control system 
and the disturbance inputs are the sum of three sinusoids, then the motion sensor outputs must 
also be the sum of three sinusoids. This is certainly true if the system is linear. Furthermore, the 
required cancellation signals must be the sum of three sinusoids in this case. Note that any sum 
of three sinusoids is specified by the six coefficients of the cosines and sines of which it is 
comprised. The ANC cancellation system is built to determine and operate on these coefficients. 

The three disturbance synchronization signals are converted in the largest analog circuit card into 
sine and cosine signals at each of the three disturbance frequencies. These are used as the 
fundamental basis functions to process both the motion sensor signals and the control commands. 

The motion sensor inputs are amplified and low-pass filtered and then are fed into the adaptive 
demodulators. The demodulators determine what combination of the three sine and cosine signal 
pairs are in each of the three motion sensor signals. For each of the three sensors there results six 
coefficients which are then fed into the A/D converters. It takes approximately 20 seconds for 
these signals to reach their final values once a new set of demodulated strut commands are issued 
by the VME computer. 

The A/D converter system, residing on one master A/D card and two slave multiplixer cards, 
converts the 18 coefficients used by the cancellation algorithm to the VME computer for 
processing. The resolution is 12 bits. 

The 12 bit D/A converters, residing on two VME cards, issue 6 commands for each of the six 
strut bending axes. The 6 commands per axis are the proportions of the 3 pairs of sine or cosine 
signals which are desired in the strut command for a single strut axis. 

The modulator multiplies the 6 coefficients per channel times the appropriate sine and cosine 
signals to compute the command for each strut’s two bending axes. 
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Motion Sensor Inputs 


Figure 6-4. Block Diagram of the ANC System 
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7. ANC Theory and Algorithms 

This chapter describes the main theoretical results that form the foundation of the ANC system. 
Section 7.1 discusses the basic problem setup and notation. The analysis of the analog 
demodulation system is described in Section 7.2. Section 7.3 covers some important basic 
results that are used to prove results in the other sections. Finally, Section 7.4 actually describes 
the algorithms used in the test visits, along with the analysis of convergence in each case. 

7.1 Control Problem Formulation and Notation 

The basic system setup is shown in Figure 7-1 below. 



rf(r) = £Re{rf°V M '} (7-D 

f=l 

iKr) = i>e{u ( 'V M '} (7-2) 

1=1 

Figure 7-1. The Basic System Notation 

The disturbance to the system and the control signal are assumed to be the sum of n sinusoids as 
shown. The input disturbance, the control vector, and the output are all assumed to be vectors of 
dimension nd, mi, and ny respectively. The system is assumed to be linear so that its output is 
also the sum of sinusoids: 

><')= |>{g(- P(- MKV*"} (7-3) 

If we make the following definitions: 

G (,) = G(-i(0> (,) (7-4) 
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III 

>■— 
e 

(7-5) 

y u) s qQ) + poyo 

(7-6) 

then we can write 



(7-7) 


The goal of the whole exercise must then be to set the expression (7-6) to zero for each of the 
three tones (z = 1,2, 3). Note that the problem for each of the three tones is entirely independent. 

So far it has been assumed that the coefficients u {,) are constants in time. The approach taken 
here, however, is to feed the computer estimates of y (,) and to allow it to change zz (,) at even 
intervals sufficiently long in duration to allow the steady state of equation (7-6) to be attained 
before the onset of a new value. For the ASTREX testbed, a duration of approximately 20 
seconds was used in most experiments. Because the control coefficients and the output 
coefficients measured at the end of the interval are time varying, a new index k must be 
introduced. 

y {i) (k) = G (i) +P {i) u (0 (k) _ (7-8) 

Also, since the problem of controlling several tones is entirely uncoupled, we drop the 
superscript (z) when considering'the problem of (7-8) for any one tone. 

7.2 Analog Demodulation System and Analysis 

We have already mentioned that the computer does not get access toy(t) directly, but only to 
estimates of the y (,> (k). We now describe how these estimates are obtained in analog electronics. 
The assumption is made that we are in the steady state so that equation (7-7) holds. Naturally, 
we must wait long enough for this to be a good approximation. 

Define: 


«.[Re{v°>} Irajv 0 ’} - 

Re{y (n) } Im{> ,(n) }] r 

(7-9) 

v(/) s [cos(a>,r) sinfcty) 

cos(co n r) sin(o) n r)j 

(7-10) 

It is then just a matter of expanding (7-7) for a single scalar element to see that 


y(t) = v(t) T q 


(7-11) 

Now construct an estimate of q as follows: 


e(t)= y(t)-v(t) T q(f) 


(7-12) 

q(t) = X v(f)e(r) 


(7-13) 


We can find out what happens to the error: 
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g(t) f A. v(t)e(t) 

A q(t) = q{t)-q 
A q(t) = -A. v(t)v(tf Aq(t) 

^\\bq(t)f = -X (v(t) T Aq(tjf < 0 


This implies that the error magnitude is monotonically nonincreasing. In fact, a little more 
analysis shows that it converges precisely to zero. Before proceeding, however, it is useful to 
show a block diagram of the demodulation system, because it is essential to a later understanding 
of the analog electronics used to implement it: 



Figure 7-2. Block Diagram of the Multi-Tone Demodulator 

Clearly, analog multipliers form an important part of this device. We now commence an analysis 
of this network from a different point of view that provides detailed information on the 
convergence behavior and sensitivity to errors. This analysis is done in the Fourier Transform 
domain with the following conventions on transforms: 

y(co)= f" = (7-15) 

27t 

It will be useful to use complex notation in what follows. To that end define: 

> ,0) (0 = ^ 2l - 1 W + jq 2i (t) (7-16) 

Keeping close track of complex operations, the estimator equations are then: 
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y(r)=£ Rc{y w (/>" M '} (7-17) 

i=! 

«(0 ® y(t) - y(0 (7-18) 

y (,) (/) = X e JW, 'e(t) + Xw i (7-19) 

where we have included an extra offset in the integrator of (7-19) because it is an important 
source of error in the electronics. We will make frequent use of the Dirac delta function in the 
following. There are several important identities required for our calculations: 

jt(f) = <? M ' <=> X (co) = 2rc8 (to - to,) (7-20a) 

x(t) = Re{v(?)} » X(( 0 ) = i(y(co) + y(-co)*) (7-20b) 

x^^y^e^ 1 <=> X (to) = y((0 - CO,) (7-20c) 


The idea is to solve for the Fourier Transforms of the coefficient estimates as a function of the 
true coefficients and the integrator offsets. To that end, we follow the signals around the loop, 
beginning with the input y(t). 

Equations (7-7) and (7-20a), (7-20b) imply: 

y«o) = (.v (,) 8 (to + co,)+ y (,r 8 (co - co,.)) (7-21) 

Equations (19) and (20c) imply: 

*(:\ 2Xtc w. . « 7i 

y ( ) (co)= — 8 (co) +—£(co-co,) (7-22) 


noting that, because e(t) is real, £(co) = £(-to)’. Equations (7-17) and (7-20b,c) imply: 

T(co) = iX(^ () ( a) + w ,)+y ( ' ) (-co-i-(o i ) ) 

Now we plug (7-22) into (7-23) and clean up the mess: 

8 (co+co,) ,8 (co-co,) 


y«o) = tai£ w ’, 


7(co+co,) ' y(co-co,) 


i-^rkco) 

“7 CO, -CO’ 


Now use (7-18) to get: 

y(co) = y(co)-£(co) = 

8 (co+co,) t 8 (co-co,) 


>-7t£ H* 


7(co+co,J 


! + £ - |£(to) = y«o) - X 7c£| H’ — 


7 (to-co,) 


£—|£(co) 

“l'CO, —CO" 1 


“7 co,- — co ■ 


8 (co+co,) ,8 (co-co,) 


7(co+co,) ' 7(co-co,) 


(7-23) 


(7-24) 


(7-24) 
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To get further define: 

ff(©).(l + £- 


Note this is zero at the tonal frequencies ± 0 ),. Consider the contribution of 

T v .Ato V 8 (©+©,) 

1 + X—- j \ w,.—7-V 

[ £?©/-CD 2 J ' _/(©+©.) 

./(©+©,)./(-©, -co,) 8 (co + co,) 

- A co, ' j (co+co,) 

- 2—5 (co + coj 

Since only the values at the zero argument of the delta function are relevant. Similarly, 

r v a® l " 1 . 8 (®-®,) 

l + \—2 - t\ w i ~~7 -r 

[ co, 2 — co 2 J y (co - co, ) 

./(co, + ©,)./(© -co,) ,8 (© — © ; ) 

A to, ' j (co — co,) 

= 2-^-8 (co-co,) 


£(co) = //(© )Y (co) - 27c ^ (w,8 (co +©,)+ w,*8 (co - co,)) 
Plugging this result into (7-22) we get: 

V (,) ((o) = — ^ H j -8 (co)+~— //(© — ©j)y(co — co ) 

J co J(D v '' 

-27t-^X( w P 5 (co-©,+©,)+w,*8 (co-co,-©,)) 

= ^//(©-©,)y(©-©,) + >.7i—8 (co-2co) 

M v '' 

_2Xtc v f.vt; . / / X\ < / / 


X - £ — 8 (©-(©, -©„))+—-—8 (©-(©,+a 

t?, ©.-©„ V ' ' p '> © + © V ' * 


Some final work is needed to compute what goes on due to the Y input. We know that H is zero 
whenever its argument is +/- an input frequency. Because Y is the sum only of such sinusoids, 
we need examine only those frequencies. The only nonzero case is for zero frequency, at which 
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-^£(co -£o,)r(co -co,)= y (,) s («>)) 
= 2rc8 (co) y (i) 


which means the coefficient is identified perfectly except for the stuff in (7-28). The final 
expression for the coefficient is: 


Y {i) (co) = 2 tt 5 (co)y (,) + X tc -^-8 (co - 2co,) 


—V —-—8 (ico — (co — co ))-i--—8 (co-(co +co 

j V ' * p V co.+co, V v 1 p. 


(7-31) 


The DC value of the coefficient is exactly right, but integrator offsets produce oscillatory 
behavior at the sum and difference frequencies of the tones, with magnitude proportional to the 
integrator gain and inversely proportional to the frequency difference. 


7.3 Helpful Theory for Cancellation Algorithms 

This section gives some essential results that are used to prove convergence for the algorithms 
that were tested on ASTREX. 

7.3.1 Standard ANC Learning Law Update Results 

In this section we consider results for the simple gradient descent algorithm with appropriate 
scaling. These results are used to show convergence for the later, more complex algorithms. 

We begin by describing the plant: 

z(k)=C(k)x, (7-32) 

where z(k) is a complex vector measurement sequence, x is a complex parameter vector to be 
estimated, and C(k) is a complex matrix. We define the following estimator: 

x(k)-. x(k + l) = x(k) + s(k) 2 C(k) H (z(k)-C(k)x{k)),x(0) = 0, 

s(k) is a real scalar such that [y(jt)|||C(ifc)|| < 1. ^ ^ 

In addition, define the following for convenience: 

A x(k) = x(k) - x 

z(k) = C(k)x{k) ( 7 . 34 ) 

e(k) = z(k) - z{k) = -C(k)Ax(k) 

With these definitions, the update equation (7-33) can be simplified, and we can write 
expressions for the parameter error update: 

x(k + 1) = x(k) + s(k) 2 C(k) H e(k) = i(k) - s(k) 2 C(k) H C(k)Ax(k) 

A x(k + 1) = Ax(k) + s(k) 2 C(k) H e(k) = (/ - s(k) 2 C(k) H C(k))Ax(k) 1 ’ 
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Before stating the main result, we start with a definition and some simple consequences. 

Definition 1 

A sequence of complex n-vectors x( k) is Li iff there exists a limit S and an M<« for every 
8 > 0 such that 0 < S - X||*(fc)|f < 8 for all r> M. In other words, ^|(x(£)|| 2 = S < «>. 

__ *=0 k=0 

Proposition 1 

The elements of an L 2 sequence converge to zero, and the sum of an L 2 sequence converges to a 
finite limit. 

Proof: 

By “converge to zero” we mean we can always choose an index in the sequenceM such that all 
elements of the sequence after it are smaller than any value 8 we care to name. 

Using the same choice of M as in the definition, 

S-thkf <8 and 0<S-XIM*)|f => lk r+1 )f <S-]£lk*)|f <8 , true for any 

*=0 k =0 *=0 

r> M. 

Also, 

|f>(t)|| S*±Wf-S<-. 

Proposition 2 

If, for a sequence of complex n-vectors X = {jc(0),x( 1),- ■ •}, Xlk^)f < a uniformly for all r, 

1=0 

then X is an L 2 sequence. 

Proof: 

The sequence of sums S r = ^||x(^)||', r = 0,1,-•• is monotonically nondecreasing and bounded 
*=o 

from above by a. Hence it converges to a bounded limit < a . 

We now proceed with the main result relating to the estimator of equation (33). These are 
deterministic results independent of conditions on the matrices C(k). 

Proposition 3 

The parameter error ||Ax(A')|| is monotonically nonincreasing, s(k)e(k ) is an L 2 sequence 
convergent to zero, and both x(£) and A x(k) converge to definite limits. 
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Proof: 

Use the scaling factor s(k) and the expression (7-35b) to obtain a bound on the new parameter 


||Ax(* +1)|| 2 = Ax(k) H (l - s(k) 2 C(k) H C(k))(I - s(k) 2 C(k) H C(k))Ax(k) 

= Ax(k) H (l - 2s(k) 2 C(k) H C(k) + s(k) 4 C(k) H C(k)C(k) H C(k))Ax(k) 

= ||Ax(A :)|| 2 - 2|j(it)C(ik)Ax(ik )|| 2 h-||^(A:) 2 C(A:)"C(A:)Ax (^)|| 2 
<||Ax(*)f -2||s(*)C(Jt)Ax(*)f +|| 5 (/:)C(/:) w || 2 || 5 (/:)C(A:)Ax(*)f 
<|Ax(*)f-2||^(Jt)C(^)Ax (/:)|| 2 H-||^(A:)C(A:)Ax (/:)|| 2 
= \\Ax(k)(-\\s(k)C(k)Ax(k)f 
Or, summarizing: 

||Ax(* + l)f ^ |Ax(*)f -|$(*)C(*)Ax(fc)f = ||Ax(A:)||' — ||^(A:)^(A:)||' (7-36) 

which cannot be an increase. The first part is proved. By iterating (7-36) we get: 

||Ax(* + of ^ ||Ax (*)|| 2 -||s(*M*)f 

N* + Of < ||Ax(* - l)f -\\s(k - 1 )e(k - l)f -|k*M*)f (7-37) 

IM + l)f < ||Ax(0 )|| 2 - Xlk"JMm)f >0 V* > 0 

m=0 

Placing the summation on the right of the last inequality establishes the second part of the 
proposition, that s(k)e(k) is an L 2 sequence and must converge to zero. 

Now use (7-35a) for the parameter estimate updates to write: 

x(k) = ^s(i) 2 C(i) H e(i) = ^(s(i)C(i) H \s(i)e(i)), 

r=0 <=o (7-38) 

||t(/)C(/) w )(5(/)e(/)|<|k/)e(0| 

By the comparison test, the parameter estimate sequence x(k ) is the sum of an Li sequence, so it 
must achieve a definite, bounded limit. So then also must the parameter error vector A x(k ), 
since it differs from the parameter estimate only by a constant. 

Notice that these results do not guarantee that the weight estimates actually converge to their 
correct values, but only that the norm of the error decreases and that the weights converge to 
some definite value. To ensure convergence to the correct weights requires that the input C(k) be 
“sufficiently rich”. The following gives one example. 
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Proposition 4 

Suppose s(k)C(k) converges to a sequence of matricesr(&) periodic with period q such that 


rank] 


n*+i) 


= nx 


U.n*+9)JJ 

Then Jc (k)->x, Ax(k)->0. 

Proof: 

Study the L 2 sequence: 

s{k)e(k) = -s(k)C(k)Ax(k) -4 -T(k)Ax„ -4 0. 

The assumption, however, guarantees that there can be no way of escaping Ax.,, = 0. 


Also, the following nonlinear equation trick can be useful in adjusting a learning stimulus. 


Proposition 5 

The equation: 

*(*) = -- bX ^ 1} .1/ . *>0, x(0)>0, *>0, a*0 

(b*+x{k-\) a Y a 

has the solution: 

x(k) = - - bX - ~- y , k> 0. 

(b a +kx(0) a ) /a 

If a > 2, then 

limXkM 2 =°° 

t = 0 

and this is not an Li sequence, though it converges to zero. 

Proof: 

First, 

bx( 0) bx( 0) 

* <0) =-- 


(b u +ox(0) a ) /a (b a y 


(7-39) 


(7-40) 
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(b a +x(k-\) a f a 


(b a +(k-\)x(Orf a 

f > 

bxj 0 ) 

(b a +(k-\)x(orf a ; 


( b * + ( k - »(or) 

(fc“+(*-l);t(0) a ) 

_ b a b a x(Q) a _ fr c ;r( 0 ) c 

’ b a (b a +(k- \)x(0) a )+b a x(0) a ~ b a + fcc( 0 ) a 
/IN bx( 0) b 


(b a +kx(0) a f a f b a 


and the proof by induction is finished for the formulae (7-39) and (7-40). Now, as to the fact that 
the sequence is not L 2 for a > 2 consider: 

^ bx( 0) (bx( 0)) _^r** 1 (bx(0)) ^ 

&\{b a +kx(0) a f a &>(b a +kx(0) a f a *=o J * (b*+kx{0) a f* 

> ±r (i,4 ° )): y dr = f (ta(0>) ’ ^ = J-_ feg): rfr 
*=° * (fc a + rx(0) a Jo (fr a + rx(0) a ^ Jo (fc a + rx(0) a f a 

(hxt mV 1 (WoiV [ , 0 S( S )C a = 2 


(^(O)) 2 r- 1 , (bx(0)} 

x( 0 ) a w s 2 / a x(0) a 


i a r else 


7.3.2 Kalman Filtering and Convergence Results 

The foundation of the Kalman Filter is in probability theory. It can be derived in numerous ways. 
Here we restate results based on Gaussian assumptions (see [ 8 ] for the proof). We use the 
brackets ( ) to denote expectation over all random variables. 
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Proposition 6 

Suppose x(k) is an nx -dimensional, complex Gaussian stochastic process with initial statistics 
(x(0)) = x 0 and & = ((*(0) - 3r 0 X*(0) - x 0 f } > 0 . Suppose w(k) and v(k) are white, Gaussian, 

zero-mean processes of dimension nx and ny, and with covariances W(k) and V(k), respectively, 
independent of each other and x. Suppose y(k) is an ny -dimensional, complex sequence of 
measurements, and that A(k) and C(k) are known complex sequences of matrices of appropriate 
dimensions to fit into the plant: 
x(k + \) = A(k)x(k) + xv(k) 

(7-41) 

y(k) = C(k)x(k) + v(k) 

Then the least mean square error estimator using the measurements is constructed as a two-step 
process, a measurement update followed by a prediction. Initialization is begun with: 


W) = * 0 
0(0) = Go 


(7-42) 


Measurement updates proceed with: 

G(*) = 0(^)-0(^)C(A:) /, (C(/:)ea)C(A:) w +V(/:))‘ , Ca)Ga) (7-43a) 


m = x(k) + Q(k)C(k) H (C(k)Q(k)C{k) H +V(k))~\y(k)-C(k)x(k)) (7 43b) 

*>0 

Predictions proceed with: 

Q(k +1) = A(k)Q(k)A(k) H +W(k) 

x(k + \) = A(k)x(k) (7-44) 

k> 0 

The variables x(k), Q(k) are the expected value of the state and its covariance given the 
measurement, and x(k), Q(k) are these just prior to measurement, respectively. 


There are a few simplifying expressions that can be used in the update equations that are useful. 
These are shown in Proposition 7. 

Proposition 7 

Suppose V(k) is invertible. Then equations (7-43) can be rewritten: 

Q(k)=(Q(ky' +c(k) H v(kr l c(k))~' (7 45 

x(k) = x(k) + Q(k)C(k) H V(ky ] (y(k)-C(k)x(k)) 
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Proof: 

Q(k)-Q(k)C(k) a (c(k)Q(k)caf+v(t)y'c(k)e(k)=(ea)-'+c(t) H v(k)-'ca)y 

because of the infamous Matrix Inversion Lemma. 

To resolve (7-43b) multiply (7-43a) on the right by C(k) H V(k)~ l to get a new expression for the 
multiplier of the measurement error: 

[e(^)-e(it)C(it) w (C(it)e(^)C(^) w +V(A:))' , C(/:)0(*)]c(/:) ,, Va')" , 

= e(it)C(Jt) w [/-(C(it)ea)C(/:) w +V(*))" , C(A:)e(^)C(/:) w ]Ka')- , 

= Q(k)C(k) H {C(k)Q(k)C(k) H +V(k))~' 

x [{C(k)Q(k)C(k) H +V(k))-C(k)Q(k)C(k)" JW 
= Q(k)C(k) H {C(k)Q(k)C{k) H +V(k)y'(V(k)y(k)- ] 

= Q(k)C(k) H (C(k)Q(k)C{k) H +V(k))'' 

This completes the proof. _ 

Next, we write some interesting results on statistics that estimate V(k). These are useful for on¬ 
line situations in which noise statistics are not known a priori. The idea is, either before or after 
updating, compute the innovation and use its square to get a handle on the noise level. 

Proposition 8 

({y(k)-C(k)x(k)Xy(k)-C(k)x(k)) H ) = C(k)Q(k)C(k)"+V(k) 

((\ik) - C(k)x(k )Xv(k) - C(k)x(k))") = V(k) - C(k)Q(k)C(k) H (7-46) 

= V(k){C(k)Q(k)C(k) H +V(k))~ l V(k) 

Proof: 

The first expression is especially easy since the noise is independent of the previous updates and 
predictions: 

(v(k){x(k)-x(k))") = 0=> 

((y(k)-C(k)x{k)\y(k)-C(k)x(k)) H ) 

= ([C(k)(x(k)-x{k))+ v(k))(C{k){x(k)-x(k))+v(k)) H ^ 

= C(k)Q(k)C(k) H +V(k) 

The second is similar but requires a little more care (use (7-45b)): 
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((j tik) - mUkf) = ((*-(*+ QC"V ~' (Cx + v- Cx))y) = -QC h => 

(( y(k ) - cU)£(fc)X;y(*) - C(*)x(*)) w ) 

= {(C(k)(x(k) - m) + v(*))((4*) - x(k)) H C(k) H + v(kf)) 

= C(k)Q(k)C(k) H - C(k)[Q(k)C(k) H ]-C(k)Q(k)C(k) H +V(k) 

= V(k)-C(k)Q(k)C(k) H 

This then is expanded (dropping the index k ) using (7-43a): 

v-cqc h = v-c{q-qc h (cqc h +v)~'cq)c h 
= V - CQC h (/ - (CQC" + v )" 1 CQC" ) 

= (cqc h +v\cqc h +v)' ] v-cqc h (cqc h +vY((cqc h +v)-cqc h ) 

= (cqc h +v\cqc h +v)~'v-cqc h (cqc h +v)'V 
=v(cqc h + v)~'v 

Now, the next job is to specialize to the case with no dynamics A(k) = 7, W(k) = 0 and the 
measurement noise is stationary at V(k) =c 2 I . The filter is restated in simplified form, without 
explicitly using the noise intensity. 

Proposition 9 

Suppose A(k) = 7, W(k) = Oand the measurement noise is stationary at V(k) = o 2 I . The filter 
can be restated as follows: 

Q(k) = (Q(k-\)-'+C(k)C(k) H y', 

cr (7-47) 

x(k) = x(k -1) + Q(k)C(k)"(y(k) - C{k)x(k -1)), Jc(-l) = * 0 

The noise estimate expressions (46) reduce to: 

((>’(*) - C(k)x(k -1 )X>’(*)" C(k)x(k -1))"} = {C(k)Q(k -1 )C{kf + iy 2 

{(ya)-C(7r)x(7:)X>’(^)-C(7:)^)) W ) = (7-C(^)j2a)C(7:) w )cf 2 (7-48) 

= {C{k)Q(k-l)C(k) H + l)~'a 2 

Proof: 

The prediction step (7-44) is trivial x{k) = x(k -1), Q(k ) = Q(k -1). The update result follows 
easily by setting Q(k) = c 2 Q(k) . The results of (7-48) follow' by these same substitutions. 
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These results are now applied to the problem of estimating a matrix that is multiplied by a time 
varying vector. For this we need the vec and kronecker product operations: 



a, i 




°ni 


a u B • 

" a \m B 

vec(A) s 

~T~ 

, A®B = 






a nX B ■ 



a nm 




vec(ABC) = (C T ® A)vec(B) 

(A®B)(D®E) = (AD®BE) (7-49) 

(A®B) r =(A T ®B t ) 

(A ® B)~' = (A ' 1 <E> B ~') 

All matrices are complex in the above, and the identities shown work only where they are well 
defined (compatible dimensions, square inverses exist, etc.). 

Proposition 10 

Consider the simplified plant of Proposition 9 except that now the plant is given by: 
y(k) = X z(k) + v(k) 

00 = 0 (7-50) 

(vec(X)vec(X) H ) = a 2 X 2 l 

with X an nyxnz complex matrix. The measurement noise has intensity V(k) = c 2 I as 
before. 

Let: 

p(k)si + Z (k) H 0(k-l)z(k) (7-51) 

Then the Kalman Filter can be expressed as: 

<!>(*) = (<D(* -1)- 1 +z{k)z(k) H Y\ O(-l) = \ 2 1 

X{k)= X(k-\) + e(k)z(k)"<t>(k), X(-1) = 0 (7-52) 

e(k) = y{k) - X(k - l)z(k) = v(k) - AX(k - l)z(k) 
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This has the solution of: 


W) = \h ' 2 / + Xz(0z(0 W ) 

X(*) = X^X^0z(0 w j+Z v (0z(0 w j^) 

AX(Jt) = X(*) - X = -X " 2 XO(^) + ^Xv(0z(0 H ^(^) 


and for v(i) s 0, 

AX(k)Q>(ky' AX(kf 

= X _4 X<t>(k)X H = A X(k- \)<t>{k - \y'AX(k -1) w 


e(k)e(k) H 

p(*) 


(7-53) 


Hence in the no noise case if X 2 » 


^£z(0z(0 H (and the ii 


inverse exists), then X(k) = X . 


This can occur afternz steps. Also, the norm of the error AX(k)<P(k)AX(k) H and its trace are 
monotonically non-decreasing. Further, e(k)/^Jp(k) is an L 2 sequence. If z(k) is bounded, 
then also e(k) = -AX (k - l)z(X) is an L 2 sequence. Finally, the rows of X ( k ) are linear 
combinations of the z(i) H , i = 0 , -,k even when noise is present. 


The noise level can be estimated using: 

((:y(*) -X(k- l)z(£)) W (y(*) -X(k -1 )z(kjj) = nyp(k)a 2 

((>’(*) - X(k)z(k)) H (yik)-X(k)z(kjf} = 


Proof: 

Rewrite (7-50a) in terms of the vec of X: 
y(k) = (z(k) T ®l)vec(X)+v(k) 


Then, identify equivalent elements with (7-4lb): 
x = vec(X) 

CU) = (zU) r ®/) 
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Expand the expressions for the filter Q: 

Q(k) = (Qtk - iy' +C(k)C(k) n y' = (d(k -1)' 1 + (z(k) T ® /„,)(;«■ ) r ® /„)") ' 

" = (fit* -1) _1 +(«(fc) r z<*)* ® 

We note that the dimension of Q is nynzxnynz, and that it begins with and maintains a 
kronecker product structure, which, on comparison with (7-52a) is just: 

Q(k) = <l>(k)'®l ny 

This is easily proven by induction if desired. Using this new expression for Q we can plug in for 
the weight update of (7-47b): 

x(k) = x(k -1) + Q(k)C(k)"(y(k)- C(k)x(k - 1)) 

= x(k -1) + ($(*)* ® /)(z(*) r ® l)"(y(k)-(z(k) T ® l)x(k - 1)) 

= x(k -1 ) + (<t>(k)' z(kY ® l)(y(k)- X(k - \)z(kj) 

U 

X(k) = X(k-l)+(y(k)-X(k-l)z(k))z(k) H ^(k) 

The formula (7-53a) is easily proven by checking its validity for k = -1 and then inductively for 
k given the result for k- 1. The second formula (7-53b) is then shown as follows: 

X(-1) = 0 

X(k ) = X(k - 1) + (y(k) - X(k - \)z(k))z(k) H <t>(k) 

= X(k- 1) + (Xz(*) + v(k) - X(k - \)z(k))z(k) H Q>(k) 

= X{k-\)(<i>(kr'-z(k)z(k) H )<l>(k) + Xz(k)z(k) H <t>(k) + v(k)z(k) H <t>(k) 

= [ x [X^(0z(0 w j + Sv(0z(0 w j^-V + Xz(0z(0 w j 

x(x ' 2/ + S^0z(0 w -za)za) w joa) + Xza0za) w W) + va)z(it) w O(it) 

= ^Xz(/)z(/) w j + Xv(0z(/) w joU) + Xza)za) w O(^) + v(^)z(^) ,/ cI)(/:) 

= zU)z(i) H j ■+ X v (*)z(0 H 
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The observations about the case when X 2 » |^z(/)z(i) w J | are obvious from the expressions 
and the size of z(i). The error expressions can be shown as follows: 


AX{k) = X(k)-X 

= x |(Sz(')z(0" ' ; / + |>0)z(0" JJo(*) + (S v(0z(0 " 

= ) + (X v(i)z(i) " j<W*) = -». ‘ ! M>(*) + (X H<)zV )" 

v(k )-0 => AX(k)Q(k)~ l AX(k)" = X~ 4 X<S>{k)X H 
= X - 4 X(<t>(k -1)" 1 +z(k)z(k) H )~'x H 

= X _4 x(<I>(jt -1) - <t>(k - l)z(k){l + z(k) H &(k - l)z(k))~'z(k) H <t>(k - l))x " 
1)X „ X- 4 X<S>(k-l)z(k)z(k) H <*(k-l)X H 


= AX(k-\mk-\y'AX(k-l) H 


AX(k-l)z(k)z(k) H AX(k-l) H 
P (k) 


= AX(k-l)<Hk-\y'AX(k-l) H 


e(k)e(k) H 
P (k) 
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The fact that the rows of X(k) are linear combinations of the z(i) H ,i = 0, -,k can be shown by 
using the rule: 

A H (al + AA H y = A H (iI-iA(l + iA‘ l A)-'iA'') 

= (l-WA(l + jA H Ay‘)iA H ={(l + iA“A)-U"A)(l + iA H Ay , U H 

= (l + iA H Ay'iA H 

A H ( a l + AA H y' =(al + A" a)~'a“ 

applied to: 

Z(*) = [z(0) - z(k)] 

V(k) = [v(0) v(*)] 

X(k) = j + X v 0')z(O w ]w) 

= [xZ(k)Z(k) H + V(k)Z(k) H \X 2 + Z(k)Z(k ) H )"' 

= [XZ(k) + V(k)]Z(k) H (X 2 +Z(k)Z(k) H )~ t 
X(k) = {(XZ(*) + V(k)\X 2 +Z(k) H Z(k))~'}z(k) H 
The statements about the convergence of the sequences are shown as follows: 
rr{AX(jt)<D(^)- 1 AX(it) w } = 

= tr{AX(* - 1 )<t>(k - 1 )" 1 AX(k - 1 )"}- — 7 7 - 
1 J P (k) 

= ,r{x4.(-l)-'X“}-iM-aO 

v*>o =*-ki£2L isL 

ti P(0 1 1 Jplkj 

Then if z(k) is bounded, 

\\z(k)\\ <z=>p(k) = \ + z(k)"<t>(k- l)z(k) <l + z 2 ||<D(A: -1)|| < 1 + r1|0(-1)|| 

Ml ^ Ml . II,,JI || AWI „ 


l + r||<D(-l)| ^|p^k)' 


>|e(*)| = ||AXa-l)za)|| is L-,. 


The noise expression (7-54a) is shown as follows, starting from (7-48): 
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tr[(y{k) -- C(k)x(k -1 jfo(k)- C(k)x(k - 1))") = tr[C(k)Q(k -1 )C(k) H + /)a 2 
(||y(Jt)- C(k)x{k - l)f ) = tr((z(k) T 0/)(0(A:)* 0 /)(z(*)‘ 0 /) + /)a 2 
= (l + z(k) T <t>(k)' z(k)' )tr(/ ny )s :2 = (l + z(k) H <t>(k)z(kj) nva 2 

This last step follows from the fact that $ is positive definite and any expression is 

therefore real (and equal to its conjugate). The other expression for the noise is derived 
analogously. 

This last proposition allows the estimation of matrices, keeping track of all algebraic structure of 
previous measurements, while using a normalized covariance matrix of order only nzxnz. This 
is a vital economy that makes the application of this algorithm for weight estimation practical, 
especially when combined with the square root filtering formulae developed in [9]. 

7.4 Adaptive Cancellation Algorithms 

This section gives the mathematical statements of the ANC algorithms that were used during 
experiments on the ASTREX testbed. The algorithms will be stated in terms of the plant 
equation (8), shown here with an additional measurement noise v(k): 

y(k) = G+ Pu(k)+v(k) 

y(k),G,v(k) e C" v , u(k)e C n ", p € C” yxnu 

7.4.1 High Pass Perturbation Algorithm 

This algorithm makes use of the fact that G is a constant, so that a high pass filter applied to (7- 
53) leaves only the Pu(k) portion of the equation, and some noise. First we define the algorithm, 
and then explore its behavior. 


Definition of High Pass Perturbation Algorithm 

Real Constants: 


A. e (0,1) 
a e (0,1] 

M max 

>° 

d(r) 


1 

low pass filter constant (-typical) 

ny(\ + nu) 

weight update constant (1 typical) 
maximum magnitude of any control component 

li *7max 

perturbation constants (<? roax = — q u = typical) 
column (r +1) of identity matrix I nu 
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Initialize: > 
r( 0) = 0 
u lo JO) = u(0) 

y lo J0)=y(0) 

u(\) = q o d(0) 

Low and High Pass Filters of y and u: 

y l0 M = y, 0 J k - D+*<(?(*)-?!■*(*-i)) 

y hi (k) = y(k)-y lm ik) 

M/oh (*) = “toJk -1) + h(y(k)-u lo Jk -1)) 
u hl (k) = y(k)-u h Jk) 

Plant Estimation: 


e P (k) = y hi (k)-P(k-\)u hi (k) 

P{k)= Ptk - 1)+ ° 

k<*)| 

e G (k) = y, m .(k)-P(k)u lm Xk)-G(k- 1) 
G(*) = G(*-l) + ae c (*) 

Perturbation Update: 




q(k) = 


r(k) = 


-q(k- 1) k even or r\(k) = 0 

- d| IksM. 

m >1 y J]{k) q > +q ( k y 


r{k - 1) k odd 
mod(r(& - 1) + l,mi) else 


Control Computation: 

u(k) = limiter(- P(k) + G^k),!*^ -q max ]+q(k)d(r(k)) 


else 


(7-56) 


(7-57) 


(7-58) 


(7-59) 


(7-60) 


The notation P(k) + means the pseudo inverse. The limiter scales its first argument such that 
each element’s magnitude is less than or equal to the second argument. 

There are several interesting things that can be said about this algorithm. The analysis will all be 
done assuming no noise ( v(k) = 0). The first step is to show that the filtering splits up the plant 
into two pieces. 
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Proposition 11 

y hi (k)=Pu hi (k) 
y lo Ak) = G+Pu lm .(k) 

Proof: 

The low pass filter operator is linear, so that 
y,oAk) = G lm (k)+Pu I0 jk ), 

G,»m = G w (*-l) + X(G-G w (*-l)), C w (0) = G 
^ Gion (k) — G =$ y, m (k) — G + Pu low (k). 

But then, by definition, 

yjk) . y(k) - y^\k) = (G+ Pu(k))-(G + Pu^(kj) 

= p(u(k) -u u Jk))= Pu hi (k). 

We now show that the High Pass Perturbation algorithm always converges to the best canceling 
control in the absence of noise. 


Proposition 12 

In the absence of measurement noise, the High Pass Perturbation algorithm produces a control 
history u(k) that converges to the optimum canceling control -P + G. Further, 

P(k) —> P, G(k) —> G. 

Proof: 

The first goal is to show that the P(k) estimator is of the form of equation (7-33) and that its 
input is “sufficiently rich” to force P(k) —» P. The first of these steps is easy, the second, very 
difficult. 


That the P(k) estimator of (7-58a,b) is in the proper form of equation (7-33) is established 
readily by setting: 

x = vec(P), x(k ) = vec(P(k)), C(k ) = u u (k) T 0 7 nv , s(k) = ■ .. 

IM*)| 


Application of Proposition 3 then ensures that P(k) converges to a fixed value, and that 
[j(fc)e(ft)||= I, ^ j |M*)| is an L 2 sequence convergent to zero. 


We now argue heuristically that the control input u hj (k ) is sufficiently rich. Equation (7-60) 
shows that the control is composed of the canceling estimate + a perturbation used to stimulate 
learning. Note that the perturbative part consists of a scaling factor^) times the periodic 
pattern d(r(k)). The pattern of the perturbation simply marches through the columns of an 
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identity matrix on every other iteration, thus ensuring that the entire space of u is spanned. 
Furthermore, the arrangement of probing first in one direction and then in the opposite direction 
on the next step with the same magnitude ensures that this attribute is carried through to u hl {k ). 
This suffices to ensure that P(k) —» P. 

Now, the G estimation system is from equations (7-58c,d) and (7-61 b): 
yioAk) = G+ Pu lo Jk) 
e c (k) = y lm (k)-P(k)u lm .(k)-G(k- 1) 

= G-G(k-\) + (P-P(k))u low (k ) —> G - G(k -1) 

G(k) = G(k - l) + a e G (k) —» G(k - l) + a {G-G(k -1)) 


The asymptotic equation for the weight G(k) ensures that it arrives at G . Then the control 
achieves its goal - P + G + the perturbation term. That the perturbation term goes to zero is seen 
from the following argument. 

First, recall that|-^^||e / ,(/:)|| converges to zero. But, the limiting of the control output ensures 
that ||m Ai (k )| is bounded by a constant, say K. However, 




M>l 


IM*>I 


ensuring that ||e / ,(^)|| itself goes to zero. Then, reviewing equation (7-59) reveals that 
asymptotically the perturbation amplitude <?(&) will obey equation (7-39) of Proposition 5 on 
every other step, so that it decays to zero. Thus the perturbation gradually disappears. 

This ends the proof. 

The conditions of equation (7-61) are reestablished after an actuator or sensor failure after several 
time constants of the low pass filter, X ~ l , have elapsed. Then the same argument again shows 
that the weights are perfectly estimated in the no noise case. 

This algorithm has the disadvantage of requiring several time constants X and a' 1 to 
converge. The next algorithm is as quick as possible, and also has better noise filtering 
capabilities. 
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7.4.2 Kalman Filter Cancellation Algorithm 

The algorithm is specified first, and then we do some analysis. 

Definition of Kalman Filter Cancellation Algorithm 
Real Constants: 

X € (0,1) low pass filter constant (—-—- typical) 

(1 + nw) 

maximum magnitude of any component of the canceling control 

M_ax 

r max » r 0 > 0 perturbation constants (r max = —, r a = 0.1 typical) 
v plant change detection factor (20 typical) 

P 2 magnitude of initial Q matrix (10 8 typical) 

d(s) column (s+1) of matrix [/„„ -I nu 0 nuxl ] 

Initialize: 

r(-l) = r 0 , 5(-l) = 2n«, d 2 (-l) = 0, ii(-\) = 2nu + \, «(-l) = 0 
e(-l) = p- 2 /, P(- 1) = 0, G(-l) = 0 
Innovation and uncertainty ratio computation: 
e(k) = y(k)~ P(k -1 )u(k)- G(k -1) 

p(k) = l + T) (k) H Q(k — l)r| (k) 


(7-62) 


(7-63) 
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Plant Change Detection and Noise Level Estimation: 
ny p (k) 

plant_changed = (a 2 (k) >v a 2 (k- 1)) a (k = 0) 
if (plant_ changed) 

et*)=(p - 1 /+n <*>n (*>")■' 


X(*) = 0 
else 


Q(k) = (Q(k-\)-'+r\(kyr\(k) H Y' 


ii = max(/7 -1,0) 


%(k) = min 1,3 


2yd 2 (k)6 2 (k -1) 
p(*)(d 2 (*)+d 2 (k - 1)) 


d 2 (k) = d 2 {k - 1) + X (d 2 (k)-6 2 (k - 1)) 


2yd \k)a 2 (k -1) 


X(*) = min 1,3 


p(£)(d 2 (k) + c 2 (k -1)) 


Plant Estimation: 

ff (k) = Q(k)i](k) 

G(k) = G{k-]) + e(k)f) l (ky (7-65 

P(k)= P(k-]) + e(k)i) 2 , nu+]) (k) H 
Perturbation Update: 

^)-(l-X(t)>_ + X») (7 . 6£ 

s{k ) = mod(s(&) + 1, 2nu+ l) 

Control Computation: 

u(k +1) = % (fc)limiter(- P(ky G(k),u nax )+ r(k)d(s(k)) (7-61 


The limiter scales its first argument such that each element’s magnitude is less than or equal to 
the second argument. 
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Before proceeding formally to the convergence results and analysis, we should cover a few points 
about the algorithm to aid in understanding its functionality. 


The variable x (k) is a heuristic measure of the degree of convergence or “confidence factor” 
in the estimation results. It is used to prevent applying control based on only a few noisy 
measurements. It varies between 0 and 1,0 representing no confidence, and 1, full 
confidence in the estimates. Before convergence it is a factor times an index showing the 
2ylc\k)o 2 (k-\) 


volatility of the noise estimate, 


(<? \k)+a\k-l)) 


, divided by the noise adjustment factor 


p (k) from equation (53) of Proposition 10. This is used in the perturbation adjustment of 
(66) and in the control computation of (67). 


• A plant change is detected if the measurement noise estimate based on the current 
measurement is very much larger than the average computed from filtering previous 
estimates (equation 64). In addition, to detect a plant change requires that a certain number 
of iterations have passed since the last plant change, as measured using the counter ii. Upon 
detecting a plant change, the normalized covariance Q(k) is set to a large value and the 
confidence factory (k) is set to zero. 


Now we turn to the main proposition concerning this algorithm, which proves that Proposition 10 
applies to this algorithm until a plant change is detected. 

Proposition 13 

Proposition 10 applies to the Kalman Filter Cancellation Algorithm with 

x=[g />], x(t)=[c(*) «*)], m=[ u ^, x-p 

until a plant change is detected. Furthermore, 

Xz(/)z0) w =Z(k)Z(k ) H , Z(k) = [z(0) - Z (k)] 

;=o 


has full rank for k > nu until a plant change is detected, no earlier than k = 2 nu + 1. 

Proof: 

It is apparent that the estimation part of the algorithm meets the conditions of Proposition 10 
with the assignments given. 

We now show that the inputs z(k), k = 0,1 , --,nu form a complete nu +1 dimensional basis. 

First, in Proposition 10, the rows of the weight estimate X(k ) are linear combinations of the 
input vectors z(i) H , i = 0,1, • • •, it. 


Thus, the P(k) have rows proportional to the u(k) H vectors input: 
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x(k) = {xz+v)(x - 2 i+z H z)-'z H 

[G(Jt) P(k)] = (XZ + V)(X- 2 I + Z H z )~' 

U0) w 


1 | M(0)" 

: i : 

i i ua)\ 


> P(k) = (XZ + V)(X- 2 1+ z H z)~ 


u(k) H 


Then the canceling control estimate - P(k) + G(k) can only have columns proportional to the 
control inputs up to that time, because the pseudo inverse preserves the spaces of its argument 
(its columns space is the row space of its argument, and visa versa). Note that the limiter in the 
control computation only scales this estimate without changing its direction. Therefore, the 
control is composed of the scaled canceling control estimate, which is a linear combination of the 
control inputs specified so far, and the perturbation term, which on examination is seen to add 
one new actuator to the scheme on each iteration: 

s(k) = k 

d{s(k )) = column k + \ of I nu 

It 


Gh[«(0) «(1) u{2) 


0 


w(m/)] = 


0 

0 


/•(0) 

0 

0 


r(l) - 
0 


Z = [z(0) z( 1) z( 2) 



0 ••• r(/iw-l) 

= upper triangular 


The * elements correspond to the canceling control estimate at each step. 

Note that the formula for r(k ) ensures that it cannot be zero for k in this range, because the value 
of x ( k ) is guaranteed to be small: 
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X(k) = min 1,3 


2y]a 2 (k)a 2 (k -1) 
p(k)(o\k)+6 2 (k-l)) 



p(k) s 1 +z(k) H <t>(k - l)z(k) = 1 + 

* ‘ 

r(k- 1) 

T* o 

OX' 2 

0 

0 


r(k - 1) 


0 

0 0 

X' 2 / 


0 


= 1 +1|*|| 2 + r (£ - 1) 2 X 2 > r(k -1) 2 X 2 


Where the * represents arbitrary values due to previous inputs. The initial r(-l) = r 0 is certainly 
not small enough to make x ( k ) appreciable. Then reference to (66) shows that the 
perturbation r(k) is near its maximum value for the times of interest. 


Thus, Z(nu ) is of full rank as required. It will remain so as long as no plant change is detected. 


The plant change is not detected unless the noise estimate (which is unbiased according to 
Proposition 10) experiences a large boost which lies several times outside of the expected 
variation of the estimate based on stationary Gaussian statistics. Also, it can not be too near a 
previous plant change. 

The “confidence factor” x (k) is fairly tolerant of variations in the estimated noise level, as long 
as they are not orders of magnitude changes, such as occur at the very start of identification, due 
to the arbitrary nature of the initial covariance choice. 

The Kalman Filter estimation algorithm exhibits the fastest possible convergence rate, because, 
as shown in Proposition 10 of Subsection 7.3.2, the exact answer is in principle obtainable in nu 
+ 1 steps of the algorithm, the number of steps required just to excite all available degrees of 
freedom in the system! The reason this is possible is because of the retention of the “covariance” 
matrix in the algorithm, which amounts to maintaining second-derivative information in the 
optimization problem, not just gradient information as occurs, for instance, in the LMS 
algorithm. 
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8. Experiments and Results 

In this chapter we discuss the various experiments run with the algorithms described in 
Chapter 7. 

8.1 ANC Trip 1 Experiments and Results 

The first ANC demonstration trip took place from Jan 8 to Jan 11, 1996. The trip was dedicated 
to the testing of the Hi-Pass Perturbation Algorithm described in Subsection 7.4.1. As mentioned 
in Chapter 2 two motion sensors were used during this test visit: linear accelerometers mounted 
in the X and Y directions on the secondary tower. Tower LPACTs 2 and 3 were used as 
disturbance sensors, and all six ACESA strut bending channels were used as cancellation 
actuators. 

Three tones were introduced at structural modal frequencies so as to cause the greatest possible 
excitation in the tower. These were at 10.85, 12.3, and 13.9 Hz. 

Open loop results are shown below in Figure 8-1. These are the power spectral densities (PSDs) 
of the voltages, after amplification, of the two accelerometers mounted on the tower. 
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It should t?e noted that the FFT measurements during this visit were made using a “boxcar” 
window on the samples. Because the signals are sines and cosines, with significant correlation 
over long times, this windowing technique caused large “skirts” around the major tones, which 
swamped out the noise at frequencies where no disturbance was present. When the dominant 
tones were canceled, this caused an apparent drop in the “noise floor,” since the skirts from the 
large tones were reduced significantly along with their parent tones. This should be kept in mind 
when reviewing the results from this first trip presented here. (Testing during the second visit 
used “Hanning” windowing, which removes this effect.) 

The algorithm was started with no prior system information and with all six cancellation actuator 
channels functioning. The algorithm was deemed to have converged when all demodulated tone 
measurements of the motion sensors were below 50 mV, the approximate noise floor at the A/D 
converters. This required about 5 minutes, and produced the results shown below in Figure 8-2. 
All closed loop results are the sum of the relevant x and y PSDs. 



-i-i-1-i-i_i_i_i_ i— v I 

10 10.5 11 11.5 12 12.5 13 13.5 14 14.5 15 

Frequency (Hz) 


Figure 8-2. Closed Loop Accelerometer PSDs for Trip 1, All Actuators On 
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The controller achieved 27 dB on the first tone, 51 dB on the second, and 26 dB on the third 
tone. This was low enough that other tone components were also present in the closed loop. 

The next step was to turn off an actuator channel, actuator 3, which corresponded to the vertical 
bending axis of strut 2 (at the 4:00 position). Convergence to a new solution required about 15 
minutes and produced the result shown below in Figure 8-3. 



Figure 8-3. Closed Loop Accelerometer PSDs for Trip 1, Actuator 3 Off 

The results were slightly improved over the case when of all actuators working. The first and 
third tones each saw about 10 dB improvement, and the second tone remained about the same. It 
is difficult to say why this was so, except that a longer period was allowed for the convergence of 
all tones, so that tones 1 and 3 converged more precisely while waiting for tone 2. At least it is 
apparent that actuator 3 was not, in itself, required for good performance. 

The next step was to turn off actuator 6, corresponding to the vertical bending axis of the 8:00 
ACES A strut. This required about 10 minutes to converge to the results shown below in 
Figure 8-4. 
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Figure 8-4. Closed Loop Accelerometer PSDs for Trip 1, Actuators 3 & 6 Off 

The results were slightly degraded from the previous case, by about 10 dB on the first tone and 
5 dB on the second. 

The most extreme failure case came next, when actuators 2 and 5 were also turned off, 
corresponding to the horizontal axis of strut 1 and the vertical axis of strut 3. Thus, only 
actuators 1 and 4 were left on, controlling the vertical axis of strut 1 and the horizontal axis of 
strut 2. Convergence was imperfect due to lack of sufficient control authority. After 30 minutes 
the result was as shown below in Figure 8-5. 
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Figure 8-5. Closed Loop Accelerometer PSDs for Trip 1, Only Actuators 1 & 4 On 

Tones 1 and 2 were still controlled down to approximately the same level as previously. Tone 3, 
however, was virtually at its open loop state. Apparently, only the third tone had insufficient 
authority to control with the available actuators. 

The results of this first trip were very good. However, there were three major improvements 
recommended: 

• Raise the disturbance levels so that visitors could see the performance improvements more 
readily. This required some changes in the amplifier gains to accommodate larger signal 
levels. 

• Create a visual, laser-based display of performance by mounting light sources on the primary 
mirror assembly and reflecting them off of the secondary and back to the primary. 

• Improve convergence speed if possible. 

These improvements were made and demonstrated on the second, final trip. 
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8.2 ANC Trip 2 Experiments and Results 

The second ANC demonstration trip took place from Jun 10 to Jun 14,1996. The trip was 
dedicated to the testing of the Kalman Filter Cancellation Algorithm described in Subsection 
7.4.2. As mentioned in Chapter 6, three motion sensors were used during this test visit: angular 
velocity sensors mounted in the X, Y, and Z directions on the secondary tower. Tower LPACTs 
2 and 3 were used as disturbance sensors, and all six ACESA strut bending channels were used 
as cancellation actuators. 

The three tonal frequencies were different from those used in the previous trip because there was 
too little output from the ACESA struts to cancel the angular motion when the input disturbance 
was turned up high enough to make a good demonstration. Transfer functions were measured 
and a new set of frequencies were chosen at 10.15, 10.625, and 12.725 Hz. 

Open loop results are shown below in Figure 8-6. These are the power spectral densities (PSDs) 
of the voltages, after amplification, of the three angular velocity sensors mounted on the tower. 



Figure 8-6. Open Loop Accelerometer PSDs for Trip 2 

The additional spurs at, say, 12.25 Hz are due to small nonlinearities in the disturbance system, 
which cause mixing of the frequencies. 
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The Kalman Filter Cancellation Algorithm was found to converge in approximately 7 minutes, as 
opposed to approximately 10 minutes for the Hi-Pass Perturbation Algorithm. The results for the 
case when all actuators are working is shown below in Figure 8-7. As before, all closed loop 
results are shown as the sum of the PSDs of the separate motion sensors. 



Figure 8-7. Closed Loop Angular Velocity PSD for Trip 2 

The algorithm achieved excellent cancellation of the tones: 33 dB at 10.15 Hz, 34 dB at 10.62 
Hz, and 44 dB at 12.72 Hz. However, the difference frequency “spurs” have been increased 
significantly, for example at 12.25 Hz. An investigation was undertaken to determine the cause 
of this effect. It was found that the commands to the actuators contained no such spurs, and that 
the sensors all showed them, including the accelerometers. Evidently, there is a nonlinearity in 
the actuation system which causes intermodulation products of the three tones to be generated. 
Even with these additional spurs, the RMS attenuation across the 9-14 Hz band was 27 dB. 

The next step was to turn off actuator 3. The new algorithm immediately detected a change in 
the plant and reconverged in 7 minutes. This is a major improvement over the Hi-Pass 
Perturbation Algorithm, which required considerably longer to converge after a failure than it did 
initially. The results for the case of actuator 3 turned off are shown below in Figure 8-8. 
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Figure 8-8. Closed Loop Angular Velocity PSD for Trip 2, Actuator 3 Off 

This case had better over all performance than the case when actuator 3 was turned on (except at 
10.15 Hz), as did the similar case during the first trip. Note also that the spur at 12.25 Hz is 
lower by about 10 dB. This would tend to suggest that Actuator 3 was the source of a significant 
portion of the nonlinearity. 

The last case is with both actuator 3 and actuator 6 turned off. Once again, convergence required 
about 7 minutes. The results are shown below in Figure 8-9. 
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Figure 8-9. Closed Loop Angular Velocity PSD for Trip 2, Actuators 3 & 6 Off 

Total attenuation was slightly worse than for the case with only actuator 3 turned off, but still 
slightly better than the full closed-loop case. Note that the spurs remain at approximately the 
same level as in the previous case, indicating that actuator 6 did not contribute much to the 
nonlinearity. 

These canceling control values were saved in a restart program as initial conditions, to allow 
quick demonstration without the need for 7 minutes of adaptation. The response of the system as 
can be judged by observing the moving laser spots was similar to turning off the disturbance 
suddenly. The system response damps out naturally over about 10-20 seconds, due to the low 
damping inherent in the structure. Of course, the cancellation system is feedforward in nature 
and does not add any damping to the open loop structure. Instead, it injects a signal that 
counteracts the effects of the disturbance at the sensor locations. 

The ANC system was adept at accommodating any actuator failure, and required no modeling 
data to operate whatsoever. 
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9. ANC Final Software Configuration 

Before discussing the construction of the software, we list the procedure for running the demo. 

9.1 Demonstration Procedure 

1. Ensure that the ASTREX testbed is floated. Ensure all actuators and sensors are powered on, 
enabled, and operational. Ensure the three tonal disturbances are turned on, and that the sync 
signals are connected to the ANC system. Ensure the lasers for visual performance 
monitoring are turned on. Turn on the power to the ANC system, and allow it to warm up for 
approx. 5 minutes. 

2. Hit carriage return several times on the user terminal, until the system is logged in and the 
OS-9 prompt appears. 

3. Type “demo” and hit carriage return. A menu of options for control are listed, as shown 
below in Table 9-1. Ensure the hardware reflects the choice chosen. If one of the restart 
options (1-3) is chosen, the program uses the saved solution to start, and then continues to 
learn. Additional actuators can be switched on or off after waiting about 1 Vi minutes. 
Reconvergence after a configuration change requires about 7 minutes. 


4. Hit control-C when you wish to exit. The actuator commands are zeroed out. The computer 
may require up to 20 seconds to respond to your control-C command. 

Table 9-1. Menu Choices for the Demo Program 


Menu Choice 

Description 

1 - All Actuators Working 

Restart from a saved solution with all actuators 
working. This gives immediate results, and then 
continues to learn. Be sure all actuators are turned on 
when this is started. 

2 - Actuator 3 Failed 

Restart from a saved solution with actuator 3 turned 
off. This gives immediate results, and then continues 
to learn. Be sure all actuators except actuator 3 are 
turned on when this is started. 

3 — Actuators 3 & 6 Failed 

Restart from a saved solution with actuators 3 and 6 
turned off. This gives immediate results, and then 
continues to learn. Be sure all actuators except 
actuators 3 and 6 are turned on when this is started. 

4 - Start Learning from Scratch 

This starts learning from scratch using any 
configuration of actuators on. It requires about 7 
minutes to converge. 

5 - Exit 

Leave the program before starting. 
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9.2 Software Hierarchy and Overview 

The ANC software is arranged in the C code structure shown below in Figure 9-1. 



Figure 9-1 Software Hierarchy 

The structure is reflected in the file names. The main routine is called “demo.c”. It handles the 
interface with the A/D and the D/A converters and the timing of the control loop. The ANC 
algorithm itself resides in “kalman2.h”. It makes use of complex matrix routines and floating 
point routines found in “floatutl.h” and “cmplx6.h”. 

We will cover this structure by first discussing the main routine, then the complex and floating 
point utilities, in preparation for the discussion of the algorithm itself. 

Throughout the software frequent use is made of data types that are pointers to structures. The 
software itself is found in Appendix A, and should be referenced in the following discussions. 

9.3 Main Routine (demo.c) 

The code for “demo.c” is arranged in distinct sections. The first section consists of the “include” 
file references. These are self-explanatory references given the hierarchy of Figure 9-1. 

The second section consists of declarations of important system constants and file names. The 
constant CHANS is the maximum number of channels that will ever be accessed for either the 
D/A or the A/D converters. Normally, the number would be 36 (6 actuators x 3 tones x 2 
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numbers (real/imaginary) per tone). However, two D/A channels had to be skipped because of 
designed-in firmware/hardware problems. 

The third code section consists of macros that are useful for manipulating or displaying data. 

The SORT and UNSORT macros map the physical channels of the hardware onto the ordered 
pairs of real numbers that the control algorithm expects. The data for this mapping is stored in 
the file “addamap.dat”. It should not be modified without due consideration, since the wiring of 
the A/D and D/A converter connectors is tied to this map. The other macros are for displaying 
the state of the control process. 

The fourth section is a little routine for handling a control-C that is typed by the user. The 
“intercept(handle_ctlC);” statement in the main routine tells the operating system to call the 
handler when the control-C is asserted. The routine sets the flag if the correct signal is asserted, 
which informs the main loop that it should quit. 

The fifth and final code section is the main routine itself. It begins with variable declarations, 
including the controllers which are of type “controller”. There is actually an array of three 
controller objects, each with its own internal state. 

Next, the A/D and D/A map is read in. Also, the table of A/D offset biases that were previously 
measured is read in. The A/D and D/A hardware is initialized. Then the controllers are 
initialized. 

Next, the user is asked what option he would like to execute. His answer is recorded as a string 
in the variable “ans”. If the answer involved using data recorded previously, the appropriate file 
is opened, and the restart data is read up using the controller routine written for that purpose. 

The control-C handler is installed. 

At last, the control loop begins. The data is retrieved and sorted from the A/D converter. The 
biases are subtracted from the readings. Then the control algorithms are called and the results are 
printed. The data is sorted for distribution to the D/A converters, and then sent. Finally, the 
controller waits the appointed time before returning to the top of the loop. 

When the control-C is detected, the loop exits and the end code is executed. First, the D/A 
converters are zeroed out. Then, the controller states are saved in the file “demo.dat”. This file 
is a restart file that could be copied over one of the files listed at the top of the code, if it 
corresponded with the correct configuration. This would be necessary if the plant were to change 
and the best cancellation results were required immediately, as for a “demo”. 

We now move on to the low level routines used for implementing complex matrix arithmetic. 

9.4 Complex Matrix Routines and LDL H Decomposition (cmplx6.h) 

Complex numbers are an integral part of the mathematics of the cancellation algorithm. 
Unfortunately, since C does not implement complex numbers as an intrinsic part of the language, 
they had to be implemented using data structures. They were implemented here using a 
structured data type: 
typedef struct 
{ float real; 
float imag; 
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} *complex; 

Thus, a complex number is a pointer to a structure with a real and an imaginary part. A complex 
number X must be allocated and deallocated, and its real and imaginary parts are accessed as X- 
>real or X->imag. 

All the usual operations are defined in the include file, whether as macros or as functions. 
Complex matrices are also defined: 

typedef struct 
{ int nrows; 
int ncols; 
complex ptr; 

} *complex_matrix; 

Clearly, a complex_matrix M must be allocated and deallocated. Its dimensions are accessed as 
M->nrows, M->ncols. A special macro is used to point to a particular complex element: 

#define cm(z,i,k) (((z)->ptr) + (i) + (k)*((z)->nrows)) 
so that cm(M,i,k) is the complex number M(i,k). Remember that indices in C start with 0! 

The final data structure defined in this include file is the LDL H decomposition. This is a 
numerically stable way of representing Hermitian positive semi-definite matrices. 

W = W h 2 0 => W=LDL h , 
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This is represented here by the following data structure: 
typedef struct 
{ int n; 
complex Lptr; 
float *Dptr; 

} *LDL_demod; 

so the L and D matrices are not represented as matrix structures, but as arrays of the right kind of 
element. If Q is an LDL decomposition, then an element (i,k) of L is accessed using a macro as 
LDL_L(Q,i,k). Similarly, the ith element of D is accessed as LDL_D(Q,i). Special routines are 
used for computing things like Q~'z to ensure that the best economies and accuracy are realized. 

9.5 Control Algorithm (kalman2.h) 

The control algorithm resides in this file. It actually is treated as an object, with its own internal 
data structure and routines for allocation, deallocation, and iteration. 
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The data structure is: 

typedef struct 
{ int ny, nu, ipert, arm; 
float lambda, sig2, r, maxu, maxr; 
complex_matrix Phat, Phatinv, Ghat, e, um, ym, eta, 

Qeta, Qetap, Qetag; 

LDL_decomp Qinv, Qinvold; 

} *controller; 

The algorithm itself is in the routine “control_iterate(y,u,c)”, where c is a variable of type 
controller. It closely follows the definition of Chapter 3, with the following variable 
identifications: 


Software 

Aleorithm 

Definition 

ym 

y(k) 

um 

u(k) 

Qinv 

Q(k )"’ 

ipert 

s(k) 

INITIAL_PERT 

r a 

(direct 

manipulation of 
perturbed actuator) 

d(s(k)) 

arm 

U(k) 

TOLER_CTL 

p- 2 


The code is intended to be commented well enough that it is relatively easy to compare with the 
algorithm definition. 

The remaining routines are just for allocation, deallocation, display, or reading from and writing 
to files. 

9.6 A/D Converter Interface (adjib.c, adjib.h) 

This is a separately compiled library for low level to access to the 601 A/D converter card, which 
also controls the A/D multiplexer boards by means of the 601 ’s channel selector. The include 
file just has prototypes for the routines that are intended to be called by users. 

The A/D board has its own on board computer. It is sent commands in ASCII via the “DPR” 
area. One sets the board up with these commands and then starts it off. The board then provides 
the data in the DPR area at the address you have told it. Experimentation was required with 
delays at certain points in the code to ensure that the system reset occurred as desired, reliably. 

Several routines just display the contents of the two status registers. 
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No interrupts are used in this code. 

This code should not be altered without close study of the manual and the existing code. 

9.7 D/A Converter Interface (dajib.c, dajib.h) 

This is a separately compiled library for low level to access to the 641 D/A converter cards. The 
include file just has prototypes for the routines that are intended to be called by users. Each D/A 
converter is separately accessed and sampled with its own routine. This was a brute force way of 
ensuring that the user did not have to know bus addresses, etc. The code is much shorter and 
more straightforward than the A/D code, because there is no on-board computer to contend with. 

Once one has studied the hardware manual a little, the code should be fairly self-explanatory. 

No interrupts are used in this code. 
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10. ANC Final Hardware Configuration 

In this chapter we briefly describe the commercial hardware that was bought for the ANC system, 
and cover in somewhat greater detail the operation of the custom analog circuit cards. For ease 
of reference. Figure 6-4 is repeated below as Figure 10-1. In the figure, the green signals are 
analog, and all others are digital. 



Motion Sensor Inputs 


Figure 10-1. Block Diagram of the ANC System 

The VME computer and peripherals, the A/D converters and multiplexers, and the D/A 
converters are all commercial off-the-shelf components. The Adaptive Demodulator and the 
Control Modulator are found on one type of custom analog card (the processing for one actuator 
and one sensor on each card), while the sensor signal conditioning and the sin/cos converter are 
mounted on the large custom analog card. Signals are passed between the analog cards on the 
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VME backplane (unused lines). The Demod Sensor Data and the Demod Strut Commands are 
passed on the D-connector terminated custom wire bundles on the front of the VME cage. 

We proceed with a discussion of the hardware from left to right in the VME card cage. 

10.1 VME Computer and Peripherals 

The VME computer is the left-most card in the card cage. It is the bus-master of the VME bus. 

It is a VME 167 single board computer manufactured by Motorola. It interfaces with the A/D and 
D/A converters via the VME bus. It interfaces with the peripherals through a small Motorola 
interface card on the back of the card cage, which plugs into the VME bus at the back behind the 
single board computer. 

Both disk drives, the floppy and the large NPI hard drive, are connected on the same SCSI bus. 
This interface is somewhat unusual for floppy drives and makes replacement of this item more 
difficult than just going to the nearest computer store. 

The NPI drive is the boot device. It has all of the software for the OS-9 real-time operating 
system (RTOS) on board. No fancy RTOS features such as semaphores, etc., are used in the 
ANC system. 

The floppy drive can be used either in the PC (FAT) format, or in the OS-9 format. The PC 
format is accessed using device /px. The OS-9 format is accessed using the /dO device. 
Remember to use the “dos” command to convert and copy files to or from /px, or else the 
carriage return/line feed sequence will not be right. 

Programming was done using the somewhat troublesome (iMACS editor. “Make” files were 
executed to compile and link, e.g. “make -f=demo.m” to recompile and link the main program 
with include files. 

10.2 A/D Converter and Multiplexers 

The A/D converter works in concert with two multiplexer boards. The A/D board is a Datel 
DVME-601B type with an ADC-12/2A module. The multiplexer boards are two Datel DVME- 
641 boards. The A/D converter is 12 bit, bipolar, with an input range of ±10 volts. 

The A/D converter is operated in the differential mode because this is the only mode which 
supports multiplexers. The multiplexer boards themselves, however, are operated in single- 
ended mode. Analog signals are only input into the multiplexer boards. Each multiplexer board 
handles up to 32 channels of input. Thus, 64 channels of input are available on the two 
multiplexer boards. Of these, only the first 36 are hooked up to the analog electronics, 32 
channels on the first multiplexer board, and 4 on the second. These appear as channels 0-35 to 
the user of the A/D interface routines. 

The 601 has an onboard Motorola MC68010 processor. Commands from the VME167 host 
computer are written into a dual-port memory area called the DPR. Data is also passed through 
this area. 

The slave multiplexer 641 boards are accessed automatically when the channel register on the 
601 board addresses a higher channel than is available on its own board. Communication 
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between the 641 and the 601 occur using the top connectors on the front of the cards, not the 
VME bus. The 641 multiplexers only take power from the VME bus. 

10.3 D/A Converters 

The D/A converters are Datel 622 analog output boards, with 12 bits resolution over a range of 
±10 volts (this is selectable). There are 16 channels per card for a total of 48 possible channels 
accessible. Of these, channel 8 (numbering from 0) on each card is unusable due to a vendor 
design problem. The controller actually requires only 36 channels. Therefore, the channel 
ranges accessed are 0-7,9-15 on each of the first two cards, and 0-5 on the third card. 

10.4 Analog Demodulation Cards 


The function and operation of the analog demodulator from an analytical point of view were 
discussed at some length in Section 7.2. At that time we presented a block diagram of the 
electronics which would be useful to reproduce here: 



Figure 10-2. Block Diagram of the Multi-Tone Demodulator 

The demodulator operates with n = 3. In addition, the cosine and sine signals are scaled, as are 
the multiply operations. We will not elaborate on the effects this has on system behavior, except 
to say that they are simple. As long as the total loop gain is maintained, the only effect is on the 
scaling of the weight estimates. 

We now describe how each of the above functions is implemented in the analog hardware. Refer 
to Appendix B for the schematics, which are on four sheets. 
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Before proceeding, it is essential to note that the multiplier is implemented with an Analog 
Devices AD734 chip, which has the behavior: 

(X\-X2)(Y\-Y2) 

10 

where all signals are in volts and the W output signal is attached to Z1 for feedback. Sheets 1, 2, 
and 3 of the schematic are nearly identical, and implement the three pairs of operations of the 
demodulator for the three distinct tones, except for the summations. For example, on Sheet 1, 
chip U1 performs the cosine multiply, U3 the sine. Chip AR1 pins 5,6, and 7 implement the 
integrator for the cosine signal, and AR1 pins 8, 9, and 10 implement the integrator for the sine 
signal. AR1 pins 1, 2, and 3 implement an amplifier to scale the q x (r) weight estimate, labeled 
as EC1 on the schematic. The final multiply prior to summing is performed by the chip U2 for 
the cosine signal, U4 for the sine signal. The sine and cosine contributions are summed in U2 
prior to being sent to the summing circuit on Sheet 4. 

The two summation operations are performed on Sheet 4, chip AR4, pins 1, 2, and 3. The sensor 
signal passes through two filtering sections prior to entering the summing junction. These are 
twin-T resonant band-pass filters with a peak gain of 3, a damping of 0.25, and center 
frequencies of 5 and 20 Hz, respectively. They produce a band pass effect between 5 and 20 Hz, 
with 40 dB/decade skirts. At the summing op amp, the gain from W5 to pin 1 of AR4 is 2, while 
the gain from any of the Y signals originating from Sheets 1-3 is 1. This relationship is due to 
the resistor values chosen in this inverting summer. The output from pin 1 of AR4 is the error 
signal E(T) on the schematic, e{t) of Figure 6-2. 

The remainder of Sheet 4 of the schematic implements the modulation of the signals emanating 
from the D/A converters. These enter as UC1 and US1 as the real and imaginary components of 
the command for the first tone, UC2 and US2 for the second tone, etc. The output of the sine 
multiply is added in as signal Z2 to the result of the cosine multiply (e.g. to U19). The result is 
fed to the inverting, gain of 1 summer of AR4, pins 8, 9, and 10. A final inverting gain change is 
used before the signal is output to the appropriate actuator (determined by soldering a jumper 
into one of W2, W3, W4, W7, W9, or W10). 

U19 is a DC/DC converter used to stabilize supply voltages on the card. 

10.5 Synchronization Signal Cards 
The synchronization signal card has two functions: 

1. Convert a pulse train into a cosine and a quadrature sine signal of the same frequency as the 
input train. 

2. Amplify and filter the sensor signals to minimize noise to the analog demodulators. Six 
channels of input are available, although only 3 are used in the ANC system. 

The schematic for this card is found in Appendix B and consists of 9 sheets. The first 4 sheets 
implement the first function. Sheets 5 and 6 implement the second function. Part of Sheet 6 and 
Sheets 7 and 8 implement processing for other sensors and are not used. Sheet 9 shows the 
DC/DC converter used to stabilize the supply voltages on the card. 
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The second function of the card, residing on Sheets 5 and 6 is discussed first because it is much 
simpler. The circuits of interest are all identical, so only the first will be discussed, 
encompassing the circuits of AR7 pins 1,2, and 3, and AR7 pins 5, 6, and 7. 

The input circuit was designed to accommodate the input from a PCB type accelerometer, which 
modulated its sensor signal on the same wires as the supply voltage, so that only two conductors 
are needed. Since that is not the function of this circuit, the diode part of the circuit was cut out, 
beginning above Cl94. The remaining circuit is a one-pole, inverting, high-pass filter with a 
comer at 1.6 Hz and a gain of 4.7 (with the value of R38 at 4.7K). The following stage, AR7 
pins 5-7, is just an inverting low pass filter with a gain of 10 and a comer frequency of approx. 

30 Hz (with Cl2 at 0.27 pF). The output of this circuit is then passed to the demodulator for the 
first sensor. The other circuits on this Sheet and the top of the next operate in an analogous 
fashion. 

The synch to cos/sin converter circuit schematics are shown on Sheets 1 through 4. It has a block 
diagram as shown in Figure 10-3 below. 


dc square 



Figure 10-3. Block Diagram of the Synch to Cos/Sin Converter Circuit 

The circuit works by creating a well formed square wave at the input frequency and filtering out 
its higher harmonics to produce a sine wave. The cosine wave is produced by shifting the phase 
of the sine wave by 90 degrees using a first order all pass. All frequency dependencies of the 
circuit are controlled to maintain phase shifts and amplitude peaks. 

The synchronization pulse train is assumed to be of TTL level but has few other known 
attributes. The first task is to create a square wave of constant amplitude and of exactly 50% 
duty cycle. This is done by using the sync signal to trigger a one-shot pulse generator, the 
duration of which is controlled by a voltage. The control voltage is developed by integrating the 
square wave output. Any deviations from 50% duty cycle create a DC bias in the average value 
of the wave. The precision limiter is used to ensure that the square wave is always switching 
between ±5 volts. Referring to Sheet 1, the variable one shot consists of the open collector 
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output, Analog Devices LM139 comparator U9, pins 1, 6, and 7, the comparator U10, and U9, 
pins 2,4, and 5 and surrounding circuitry. The bias integrator is AR1, pins 12, 13, and 14. The 
precision limiter is constructed from the analog switches of U3. 

The next stage of the circuit is the voltage-controlled band-pass filter, the output of which is the 
sine wave, on pin 14 of AR2. This is essentially a twin-T two-pole resonant band-pass, the 
frequency of which is controlled with what amounts to a voltage-controlled resistor formed by 
the analog multiplier U12. The control voltage emanates from the integrator output, pin 7 of 
AR2. The square wave and the cosine are nominally in quadrature, so that their product should 
produce no DC bias. However, phase shift errors in the sine wave resulting from improper 
centering of the band-pass filter frequency response must produce similar shifts in the cosine 
wave, resulting in a DC bias for the product of the cosine and the square waves. This product is 
created by U13 at pin 12, integrated, and used to correct the center frequency of the band-pass. 


Finally, on Sheet 4, the sine wave is passed through a one-pole all pass filter 


(*-©„) 


(5 + 0 ) „) 

produce the cosine wave. This filter has a gain of 1 over all frequency, but its phase varies from - 
180 degrees at DC to 0 degrees at infinity. It attains -90 degrees at the frequency s = j(£> 0 . This 
is the desired situation. The center frequency is adjusted with another voltage controlled resistor 
consisting of U25, with an input coming from the integrator output at pin 1 of AR2. The 
integrator has as its input the product of the sine and cosine waves, which has zero DC bias if the 
phase between the sine and cosine is maintained at 90 degrees. Variations in the phase cause a 
DC bias which integrates to shift the center frequency of the all-pass until the phase is corrected. 


Two more circuits identical with the one just described are implemented for the remaining two 
tones on Sheets 2, 3, and 4. 


This clever circuit was designed by Mr. Bill Burton, a principal electrical engineer at Harris 
GASD. 
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10.6 Custom Connector Wire List 

The custom D-connectors that connect the analog demodulator cards with the A/D and D/A 
converters have a wire list as shown below in Table 10-1. The signals are grouped by 
demodulator card. Each demodulator card has one demodulator for three tones, corresponding to 
one sensor, and one modulator for three tones, corresponding to one actuator. The signals are 
grouped by demodulator. The terminology “y cos 1” indicates that this is the cosine component 
of tone 1 for this demodulator card’s sensor measurement. Similarly, “u sin 2” indicates that this 
is the sine component of tone 2 for this demodulator card’s actuator command. 

The pin number on the demodulator card’s D connector is the second column. The destination 
card for the wire is listed in column 3. The label on the corresponding connector is in column 4. 
The pin number on the connector is in column 5, and D/A channel or A/D channel of interest is 
listed in the remaining two columns. 

Table 10-1. Wire List for the Custom Connector 


Signal 

Demod 

Pin 

Dest. 

Card 

Connect. 

No. 

Connect. 

Pin 

D/A 

Channel 

A/D 

Channel 

Demod 1 







y cos 1 

1 

Mux 1 

P3 

24 


0 

y sin 1 

2 

Mux 1 

P3 

10 


1 

y cos 2 

3 

Mux 1 

P3 

21 


2 

y sin 2 

4 

Mux 1 

P3 

7 


3 1 

y cos 3 

5 

Mux 1 

P3 

18 


4 

y sin 3 

6 

Mux 1 

P3 

4 


5 

u cos 1 

7 

D/A 1 

P9 

24 

0 


u sin 1 

8 

D/A 1 

P9 

10 

1 


u cos 2 

9 

D/A 1 

P9 

21 

2 


u sin 2 

10 

D/A 1 

P9 

7 

3 


u cos 3 

11 

D/A 1 

P9 

18 

4 


u sin 3 

12 

D/A 1 

P9 

4 

5 


a. return 

13 






Demod 2 

Demod 

Pin 

Dest. 

Card 

Connect. 

No. 

Connect. 

Pin 

D/A 

Channel 

A/D 

Channel 

y cos 1 

1 

Mux 1 

P3 

15 


6 

y sin 1 

2 

Mux 1 

P3 

1 


7 


3 

Mux 1 

P2 

24 


8 
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Demod 6 

Demod 

Pin 

Dest. 

Card 

Connect. 

No. 

Connect. 

Pin 

D/A 

Channel 

A/D 

Channel 

y cos 1 

1 

Mux 1 

P3 

3 


30 

y sin 1 

2 

Mux 1 

P3 

14 


31 

y cos 2 

3 

Mux 2 

P4 

24 


32 

y sin 2 

4 

Mux 2 

P4 

10 


33 

y cos 3 

5 

Mux 2 

P4 

21 


34 

y sin 3 

6 

Mux 2 

P4 

7 


35 

u cos 1 

7 

D/A 3 

Pll 

24 

0 


u sin 1 

8 

D/A 3 

Pll 

10 

1 


u cos 2 

9 

D/A 3 

Pll 

21 

2 


u sin 2 

10 

D/A 3 

Pll 

7 

3 


u cos 3 

11 

D/A 3 

Pll 

18 

4 


u sin 3 

12 

D/A 3 

Pll 

4 

5 


a. return 

13 







10.7 Bias Adjustment Procedures 

It is quite likely that the integrator biases will need to be adjusted out on the demodulator cards 
from time to time. In addition, the bias table used by the software may need to be updated 
occasionally, certainly every time the bias adjustment potentiometers are zeroed, or when the bias 
table measurements of Subsection 10.8.2 show biases greater than about 50 mV. The following 
gives the two procedures to use in these cases. 

10.7.1 Adjusting Out the Integrator Bias on the Demodulator Cards 

This procedure should be followed when the biases measured in the procedure of 10.8.2 get 
larger than about 50 mV. 

1. Turn the ANC system power on. Ensure valid sync signals are being input for all three tones. 
Ensure that sensor signals are disconnected or turned off. Ensure that the actuators are 
disabled or turned off. Wait 5 minutes for the electronics to stabilize. Ensure all the 
switches on the demodulator cards (see Figure 10-4 below) are in the “up” position, as they 
should be for operation. 

2. Run the program “chkdaad” by typing this as a command at the OS-9 prompt. Respond to 
prompts for files. The first file required is the map for the A/D and D/A conveners, 
“addamap.dat”. The second file is the data file to use for commanding the D/A conveners, 
“chkdaad.dat”. The program prints out the measured A/D readings every 5 seconds. These 
should be zero, but are normally small values less than 50 mV. The sin and cosine sensor 
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readings for the various tones are shown with the corresponding channel number displayed by 
the chkdaad computer program below in Table 10-2. 

3. Using a small screwdriver, adjust each of the potentiometers corresponding to channels with 
large biases to less than 20 mV or so. The potentiometers as viewed from the front of the 
VME cage are shown below in Figure 10-4. Counterclockwise rotation of the potentiometers 
makes the biases more positive. Note that the tone 1 sine potentiometer is above the tone 1 
cosine potentiometer, in contrast with the order for the other tones shown in Figure 10-4. Hit 
control-C to exit the chkdaad program when finished. 

4. Rerun the bias table computation of Subsection 10.7.2. 


Table 10-2. Channel Map for the CHKDAAD Program 


Demod 1 

A/D Channel 

Demod 4 

A/D Channel 

tone 1 cosine 

0 

tone 1 cosine 

6 

tone 1 sine 

1 

tone 1 sine 

7 

tone 2 cosine 

12 

tone 2 cosine 

18 

tone 2 sine 

13 

tone 2 sine 

19 

tone 3 cosine 

24 

tone 3 cosine 

30 

tone 3 sine 

25 

tone 3 sine 

31 

Demod 2 

A/D Channel 

Demod 5 

A/D Channel 

tone 1 cosine 

2 

tone 1 cosine 

8 

tone 1 sine 

3 

tone 1 sine 

9 

tone 2 cosine 

14 

tone 2 cosine 

20 

tone 2 sine 

15 

tone 2 sine 

21 

tone 3 cosine 

26 

tone 3 cosine 

32 

tone 3 sine 

27 

tone 3 sine 

33 

Demod 3 

A/D Channel 

Demod 6 

A/D Channel 

tone 1 cosine 

4 

tone 1 cosine 

10 

tone 1 sine 

5 

tone 1 sine 

11 

tone 2 cosine 

16 

tone 2 cosine 

22 

tone 2 sine 

17 

tone 2 sine 

23 

tone 3 cosine 

28 

tone 3 cosine 

34 

tone 3 sine 

29 

tone 3 sine 

35 
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tone 1 sine 
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tone 1 cosine 
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tone 2 cosine 
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tone 2 sine 
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tone 3 cosine 
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®l 

tone 3 sine 


Figure 10-4. Arrangement of the Potentiometers on the Demod. Cards 
10.7.2 Recomputing the Bias Table 

This procedure should be run at least every time the potentiometers are readjusted on the 
demodulator cards. 

1. Turn the ANC system power on. Ensure valid sync signals are being input for all three tones. 
Ensure that sensor signals are disconnected or turned off. Ensure that the actuators are 
disabled or turned off. Wait 5 minutes for the electronics to stabilize. Ensure all the 
switches on the demodulator cards (see Figure 6-4 above) are in the “up” position, as they 
should be for operation. 

2. Run the program “measbias” by typing this as a command at the OS-9 prompt. Let the 
program run for about 5 minutes to get good average values for the biases. Hit control-C to 
exit. The new bias table readings are stored in the file “measbias.dat,” where they are read by 
the “demo” program. 
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11. Conclusions and Recommendations for Future Work 

This report has described efforts undertaken under the Adaptive Neural Control (ANC) Program. 
Previous developments had theoretically and experimentally demonstrated neural adaptive 
feedback control and feedforward control for tonal disturbances. The ANC effort rounds out this 
development by first demonstrating feedforward control for the difficult case of broadband, 
continuous spectrum disturbances and then integrating and demonstrating a multi-tone, fault- 
tolerant cancellation in the ASTREX testbed. Thus the demonstrated capabilities cover all 
disturbance types likely to be of relevance in applications. In addition to the adaptive control of 
flexible space systems which was the focus of this effort the ANC architecture lends itself to 
numerous dual-use and commercial applications. In the remainder of this chapter we describe 
some of these dual-use applications and make recommendations for broader efforts to advance 
existing ANC capabilities to support these applications. 

A key feature of the ANC architecture is its hierarchy of standardized components. This 
modularity permits great implementation flexibility and also facilitates a multiplicity of diverse 
applications. As Figure 11-1 illustrates, only a relatively small number of basic components need 
to be developed and perfected and then combined to support an enormous number of 
applications. This feature greatly reduces nonrecurring engineering costs and broadens the scope 
of dual-use activities. 

In the area of noise and vibration control, potential applications of ANC range from vibration 
isolation of precision space sensors to consumer products. In the body of the report, we have 
already alluded to the use of ANC for autonomous adaptive vibration isolation of space payloads 
from spacecraft—generated disturbances. There are also many industrial applications including 
the quieting of machinery in manufacturing facilities and the isolation of semiconductor 
production facilities from ground vibration. 

In addition, we have identified numerous applications for relatively low cost, high volume 
consumer products. Table 11-1 lists some of the possibilities. As well as noise and vibration 
applications involving the use of both our commercial actuator technologies and ANC, there are 
systems that would employ only ANC adaptive controllers. These items include heating and air- 
conditioning regulators, appliance power regulators and active automobile suspension controls 
that use a neural controller for fault identification and controller recovery. 

Of even broader scope are the applications of ANC outside of noise and vibration control. Some 
of these are listed in Table 11-2. These items involve ANC in both an adaptive control role or as 
a passive damage or fault detection system. The latter application would be useful, for example, 
in the detection and location of damage in civil engineering structures, especially in areas of high 
seismic activity. Of particular note are two specifically identified applications to medical 
diagnosis; using passive acoustical measurements to detect failures in Bjork Shiley coronary 
values and to diagnose coronary occlusions. Although, in the former case, some work has been 
done using a simple multilayer perceptron to give a binary output decision, the general ANC 
architecture has the potential for greatly improved diagnostic accuracy. Because it can 
incorporate detailed models of dynamic processes, the ANC architecture offers new dimensions 
that go well beyond "black box" neural classifiers. These new capabilities could include the 
ability to fold in directionality and acoustic propagation dynamics information using multiple 
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sensors, the incorporation of dynamic models of valve operation and mechanical vibration and 
the use of these internal models to filter out variabilities in the data that are irrelevant. In other 
words, an ANC system can accept and then refine internal dynamical models to progressively 
sharpen its diagnostic accuracy. 

An "Ecosystem" of Intelligent System Products 


Intelligent 

Mechanical 

Products 


Intermediate^ 

Components 


Manufacturing 

Active noise 

Aerodynamic 

robots 

cancellation 

controls 

Robotic vehicles 

Vibration 

Autonomous 

Robot planetary 

control 

spacecraft 

explorers 

Design 

Management 


optimization 

planning 

Adaptive servo 
positioners 

systems 

systems 



Figure 11-1: Completely Modular, Standardized ANC Components Support a Broad- 
Based Industry - From Basic Module Supply to Individual Applications 
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Table 11-1: Potential Consumer Applications of Harris’ Vibration Control Hardware and 
Adaptive Neural Algorithms 

Actuator Technology - provides compact, unobtrusive actuation with maximum cost- 
effectiveness 

Neural Adaptive Control - provides "intelligent" control for noise and vibration reduction. 

Able to work autonomously and adapt to rapidly changing 
conditions. 


Small Actuators + Neural Control 


System 

Benefits 

Washing Machines 

less noise, less wear, longer life 

Dishwashers 

less noise, less wear 

Quieter Lawn Mowers, 

reduces engine vibration (acoustic transducers and neural 

Weed Trimmers, etc. 

controls reduce far field emissions) 


Larger Actuators + Neural Control 


System 

Benefits 

Automobile Engine Mounts 

eliminates transmission of engine vibration into 
passenger cabin 

Outboard Motors 

reduces acoustic noise, hull vibration 


Neural Control with Acoustic Transducers 


System 

Benefits 

Vacuum Cleaners 

obnoxious noise eliminated for user and bystanders 

Hair Dryers 

reduces distant noise, reduces interference with other 
activities (conversation, TV. etc.) 

Quieter Lawn Mowers, Weed Trimmers, 

reduces far field acoustic emissions, improves user 

etc. 

comfort and neighborhood tranquillity 


Neural Controls 


System 

Benefits 

Heating/Air-Conditioning 

maintains desired temperature while minimizing power 
consumption 

Appliance power regulators (outlet 
attachment, neural control and appliance 
transducers regulate power consumption) 

reduced electricity bills 

Active automobile suspension control 
(neural control for fault ID and controller 
recovery) 

system fault tolerance and reliability 
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Table 11-2: Adaptive Neural Control Applications Outside Noise & Vibration Control 
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Of course, most of the commercial and medical applications sketched above are well into the 
future. A great deal needs to be done to perfect ANC systems, even for the space or airborne 
DoD applications. 

One important area of development is implementation refinement. To date for the most part, we 
have been using existing DSP hardware. This is appropriate for the early stages of development 
when it is necessary to accumulate sufficient experience in the operational characteristics of 
ANC algorithms to be able to formulate meaningful performance and design specifications for 
improved, dedicated hardware. We are nearly at the point at which such specifications can be 
elaborated. It is clear that there are two distinct directions that can be taken in hardware 
improvement. The first is to develop analog implementations. Although the simplest algorithms 
have been demonstrated with analog hardware, this remains to be done for the more complex 
systems. Analog implementation eliminates the need for a computer and associated software and 
reduces cost for many applications. Since the more sophisticated ANC systems will still require 
digital capability, the second direction for advancement is the development of ASICs for specific 
applications. Such customized ICs would likely embody a number of standardized, higher-level 
modules in the ANC architecture. Properly designed, these could serve as the basic building 
blocks from which a variety of diverse systems could be built. 

A second broad area for development deals with the refined integration of ANC hardware with 
actuator hardware and with structural components. Harris has, over the past decade, developed 
several lines of actuator hardware (both inertial and intrastructural). A natural extension is to 
integrate the intelligent control capability of ANC directly into each (collocated) actuator/sensor 
unit. One obtains structures that are "smart" because the actuators, joints and fittings are 
intelligent. This "smart actuator" concept relies on the fact that even the global, centralized ANC 
model reference adaptive control algorithm can be re-partitioned into parallel components, each 
of which resides within one collocated sensor-actuator unit. Figure 11-2 shows the basic 
concept. The forward and backward path signal flows associated with the particular sensor and 
actuator unit are handled locally. At the same time, weight adjustment computations associated 
with interactions among the sensor/actuator units are executed using forward and backward 
signal information from the other sensor/actuator units transmitted through a data bus. Thus, all 
computations are carried out within the local processors, yet by virtue of the inter-processor 
communication, all processors (that are in communication) work in combination to address the 
multi-input, multi-output centralized control problem. 

A further possibility is the integration of ANC technology into smart structures involving 
embedded "mesoscale" actuation and sensing devices dispersed within the structural material. 
Clearly, the development of mesoscale actuators realizes the desirable feature of highly 
redundant and distributed sensory and actuation devices. However, to properly exploit this 
technology one needs some form of controller that is adaptable, autonomous, highly parallel and 
fault tolerant. For this reason, the implementation of massively parallel, decentralized adaptive 
control algorithms (in forms that can be explained in terms of neural network algorithms) seems 
a perfect complement to mesoscale actuation technology. 


132 


ANC Final Report 


Intelligent Actuation 

- processing is distributed among the actuator/sensor units 

- all the actuator/sensor units that are in communication 
with one another work together to solve (and implement) 
the optimal centralized control law 



One unit: converges 
to the optimal SISO 
control law 


Two or more 
noncommunicating 
units converge to the 
optimal decentralized 
control (with local 
SISO loops) 


Two or more units 
communicating 
through a data bus 
converge to the 
optimal centralized. 
MIMO control 


But, as in all previous 
cases, all processing 
still occurs within the 
individual units 


Figure 11-2: SMARTACT’s (Smart Actuators) 
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The ANC. systems considered to-date all involve a centralized information pattern—i.e. full 
connectivity among all the sensors and actuators. However, in systems having a large number of 
small actuation and sensing devices spatially distributed throughout the structural matrix, 
centralized control entails an impractically large number of interconnections. Thus, for 
application of adaptive neural control to mesoscale actuation systems some form of 
decentralization constraint needs to be imposed on the controller structure. There are a variety of 
ANC structures for implementing this kind of decentralized control, with suitable limited 
interconnection paths, that is needed for coordination of a large number of sensors and actuators. 

The most straightforward approach involves the basic ANC controllers described above, but with 
interconnection constraints. What this means is that some of the synaptic connections that would 
otherwise exist in a centralized system are removed, leaving only a subset. The remaining 
synaptic connections correspond, one-to-one, to a decentralized information pattern, e.g. each 
actuator being driven by a subset of sensors that are spatially proximate to the actuator. 

This pattern of interconnections is illustrated by Figure 11-3. The decentralized learning inherent 
in ANC readily accommodates such decentralization in the control-indeed none of the 
fundamental stability and convergence properties are modified. Decentralized information 
patterns simply mean that only a subset of synaptic connections are active and with respect to the 
corresponding weights, the system is still a gradient descent machine. A decentralized ANC 
controller still seeks a least mean-square error condition, but operates under the constraints 
imposed by decentralized structure. When a decentralized control structure is imposed, then the 
controller can obviously be decomposed into independent sub-controllers-for example, one such 
subcontroller would activate one actuator, coordinating its set of geographically proximate 
sensors. This would permit the use of small, computationally simple processing units that are 
spatially distributed throughout the structure. 
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• Embedded sensors and actuators widely dispersed within the structural 
material 

• Co-location not necessary 

• Parallel neural processors - each one associated with one actuator - each 
processor very simple (few neurons) 

• Each actuator processor connected with a spatially local subset of sensors 

• Without detailed model information, processors implement adaptive 
vibration damping and disturbance cancellation 

Figure 11-3: Distributed Adaptive Neural Architecture for Intelligent Structures 
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Appendix A. C Code 


Main Routine (demo.c) 


/* Main ANC control task. Kalman filter algorithm with state save. */ 


#include <stdio.h> 

#include <math.h> 

#include "ad_lib.h" 

#include "da_lib.h" 

#include "floatutl.h" 

#include "cmplxS.h" 

#include n kalman2.h” 

/*====================== IMPORTANT CONSTANTS =====================*/ 

#define CHANS (38) 

#define SAMPLE_PERIOD (20.0) 

#define ULIMIT (5.0) 

#define NU (6) 

#define NY (3) 

#define MAXN (6) 

char *filenames[] = {"all_ok.dat", 

"fail_3.dat", 

-fail_3_6.dat"); 

/*=========================== MACROS ===============================*/ 

#define trunc(v) (l.e-3 * floor((v)*1.e3)) 

#define SORT(dst # src.vec,cnt) \ 

{ for (i = 0; i < cnt; i++) \ 

{ dst[i] = src[vec[i]]; \ 

> \ 

> 

#define UNSORT(dst,src,vec,cnt) \ 

{ for (i = 0; i < cnt; i++) \ 

{ dst[vec[i]] = src[i]; \ 

} \ 

} 

#define PRINT_STEP_HEADER \ 

printf("\n====== New Step, actuators %d to %d, sensors %d to %d 


======", \ 

0, (NU-1), 0, (NY-1)); 

#define PRINT_Y_U \ 

{ printf (" \n-controller %d: - 

- n ,k); \ 

printf("\ny: "); \ 

for (i = 0, 1 = 0, j = 2*k*MAXN; i < NY; i++, j+=2) \ 


{ if (1 < 4) { 1++; } else (1=0; printf("\n "); } \ 
printf(" (%7.3f ,%1 .3f)",trunc(yytj]),trunc(yy[j+1])); \ 


/* seconds */ 

/* volts maximum allowed out */ 
/* number of control outputs */ 
/* number of measurements */ 
/* max NU or NY */ 
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} \ 

print'f (" \nu: "); \ 

for (i - 0, 1 ■ 0, j - 2*k*MAXN; i < NU; i++, j+=2) \ 

{ if (1 < 3 ) { 1 ++; } else { 1 = 0; printf("\n "); ) \ 
printf(" (V7.3f,%7.3f)",trunc(uu[ 3 ]),trunc(uu[j+1])); \ 

} \ 

if (controllers[k]->ipert < NU) \ 

{ printf("\nperturbation of %f at control input %d", \ 
controllers[k]->r,controllers[k]->ipert); \ 

} \ 

else if (controllers[k]->ipert < 2*NU) \ 

{ printf("\nperturbation of %i at control input %d", \ 
-controllers[k]->r # controllers[k]->ipert-NU); \ 

) \ 

else \ 

{ printf("\nno perturbation"); \ 

) \ 

} 

/*sssssbsss==ssbbssbssbb Control C Handler ==================*/ 

int ctlC_flag «= 0; 

int handle_ct 1 C(signal) 
int signal; 

{ ctlC_flag = (signal < 32); 

) 

/*====bss==================== MAIN ==========================*/ 

main() 

{ 

int i, 3 , k # 1 , sleep_count, ok; 

FILE *fp; 
char ans[5]; 

int advec[CHANS], davec[CHANS]; 
float adbuff[CHANS], dabuff[CHANS] # uu[CHANS], yy[CHANS]; 
float yybias[CHANS]; 

controller controllers[3]; 

/* Convert from seconds to sleep count */ 

sleep_count = (int) (SAMPLE_PERIOD * 100.0); 

/************************************************************/ 
/* initialization of da and ad mapping vectors advec, davec */ 
/a***********************************************************/ 

fp b fopen("addamap.dat","r"); 

for (i = 0; i < CHANS; i++) 

{ fscanf (fp,"S 6 d" / tadvec[i]); 

} 

for (i = 0; i < CHANS; i++) 

{ fscanf(fp,"%d" / &davec[i]); 
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dabuffli] = 0.0; uu[i] = 0.0; 

> •' 

fclose(fp); 

fp = fopen("measbias.dat","r"); 

for (i = 0; i < CHANS; i++) 

{ fscanf(fp,"%e",&yybias[i]); 

> 

fclose(fp); 

/a*************************************/ 

/* A/D and D/A Initialization */ 

/♦♦a***********************************/ 

ad_setup(0,CHANS); 
da0_setup(); 
dal_setup(); 
da2_setup(); 

/**************************************/ 

/* Controller Initialization */ 

controllers[0] = control_new(NY,NU,ULIMIT); 
controllers[1] = control_new(NY,NU # ULIMIT); 
controllers[2] = control_new(NY,NU,ULIMIT); 

/* Load data for restart. */ 

do 

{ printf("\n"); 

printf("\nRestart Option (1 to 3):"); 
printf("\nl -- All Actuators Working"); 
printf("\n2 — Actuator 3 Failed"); 
printf("\n3 -- Actuators 3 & 6 Failed"); 
printf("\n4 -- Start Learning from Scratch"); 
printf("\n5 — Exit \n> "); 

gets(ans,4); 

> while (ans[0] < '!• || ans[0] > '5'); 

if (ans[0] >= '1' && ans[0] <= '3') 

{ fp = fopen(filenames[*ans- , l'],"r"); 

control_restart(fp,controllers[0]); 
control_restart(fp,controllers[1]); 
control_restart(fp,controllers[2]); 

fclose(fp); 

printf("\n\nRestarting from stored data in file <s&sx\n\n" 
filenames[ans[0] - 'l']); 

> 

else if (ans[0] == '5') 
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{ printf("\n"); 
exit(0); 

} 

/* install handler for A C */ 
intercept(handle_ctlC); 

/♦•a***********************************/ 
/* Main control loop */ 

/**************************************/ 

while(!ctlC_flag) 

{ 

/* Get A/D converter data */ 
ad_convert{adbuff,CHANS); 


/* Sort adbuff into yy. Apply control algorithm(s). */ 

/* Unsort uu bach into dabuff. */ 

/*=====s:s:=s:=== = = = = = = = = = = = = = = = = = = = = = = = ssssbsbsbsssssssss*/ 


SORT(yy,adbuff,advec , CHANS); 

for (i = 0; i < CHANS; i++) 

{ yy[i] -= yybias[i]; 

> 

PRINT_STEP_HEADER; 

control_iterate(tyy[0] , tuu[0] , controllers[0]); 
k = 0; PRINT_Y_U; 

control_iterate(&yy[12] , tuu[12] # controllers[1]); 
k = 1; PRINT_Y_U; 

control_iterate(&yy[24] ,&uu[24],controllers[2]); 
k = 2; PRINT_Y_U; 

printf("\n"); 

UNSORT(dabuf f,uu,davec,CHANS); 

/*******************************/ 

/* Output D/A Converter Data */ 
/*******************************/ 


da0_convert(dabuff,0,16); 
dal_convert(fidabuff[16],0,16); 
da2_convert(fitdabuff[32] ,0,(CHANS-32)); 

/* Wait the requisite time before taking more data. */ 

tsleep(sleep_count); 


/* Zero out the D/A converters. */ 
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for (i = 0; i < CHANS; i++) 
dabuffti] = 0.0; 

daO_convert(dabuff,0,16); 
dal_convert(fcdabuf f[16],0,16); 
da2_convert(fcdabuff[32],0,(CHANS-32)); 

printf("\nda converters have been zeroed! \n\n"); 

/* Save data for later restart. */ 

fp = fopen("demo.dat","w"); 
control_save(fp,controllers[0]); 
control_save(fp,controllers[1]); 
control_save(fp,controllers[2]); 
fclose(fp); 

} 
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Kalman Filter Cancellation Algorithm (kalmanl.h) 

I **************************************************** I 

/* This file implements the HIMO single tone */ 

/* algorithm based on optimal filtering of y. */ 

/****************************************************/ 

#undef FALSE 
#define FALSE (0) 

#undef TRUE 
#define TRUE (1) 

#define TOLER_CTL ((float) 1.0e-8) 

#define INITIAL_PERT ((float) 0.1) 

#define MIN_PN0RM ((float) 0.1) 

#define ARM_COUNT(nu) (2*(nu+l)) 

/************************************************/ 

/* controller type definition */ 

/************************************************/ 

typedef struct 
{ int ny, nu, ipert, arm; 

float lambda, sig2, r, maxu, maxr; 

complex_matrix Phat, Phatinv, Ghat, e, urn, ym, eta, 
Qeta, Qetap, Qetag; 

LDL_decomp Qinv, Qinvold; 

} *controller; 

/************************************************/ 

/* New controller allocation */ 

/************************************************/ 
controller control_new(ny,nu,ulimit) 
int ny, nu; 
float ulimit; 

/****************** ******************************/ 

{ controller c; 

c - (controller) malloc(sizeof(*c)); 

c->nu = nu; 
c->ny = ny; 
c->ipert = 2*nu; 
c->arm = ARM_COUNT(nu); 

c->lambda = 1.0 / ((float) (1 + nu)); 
c->sig2 = 0.0; 

C->r = INITIAL_PERT; 
c->maxu =0.9 * ulimit; 
c->maxr = 0.1 * ulimit; 

c->Phat = cm_new(ny,nu,0.0,0.0); 
c->Phatinv = cm_new(nu,ny,0.0,0.0); 
c->Ghat = cm_new(ny,1,0.0,0.0); 
c->e = cm_new(ny,1,0.0,0.0); 
c->um = cm_new(nu,1,0.0,0.0); 
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c->ym = cm_new(ny,1,0.0,0.0); 
c->eta = cm_new(nu+l,1,0.0,0.0); 
c->Qeta = can_new(nu+l, 1, 0.0,0.0) ; 
c->Qetap = cm_new(nu,1,0.0,0.0); 
c->Qetag = cm_new(l,1,0.0,0.0); 

c->Qinv = LDL_new(nu+l); 
c->Qinvold = LDL_new(nu+l); 

LDL_eye(c->Qinv,TOLER_CTL); 

return c; 

} /* control_new */ 

/A***********************************/ 

/* Controller Deallocation */ 

/************************************/ 
void control_£ree(c) 
controller c; 

/************************************/ 

{ LDL_£ree(c->Qinvold); 

LDL_free(c->Qinv); 

cm_free(c->Qetag); 
cm_free(c->Qetap); 
cm_free(c->Qeta); 
cm_free(c->eta); 
cmfree(c->ym); 
cm free(c->um); 
cm_free(c->e); 
cm_free(c->Ghat); 
cm_free(c->Phatinv); 
cm_free(c->Phat); 

free(c); 

) /* end control_free */ 

/************************************************/ 

/* Controller restart */ 

/* c is already allocated and zeroed. */ 

/* fp is open for read, pointing at ghat */ 

/************************************************/ 
void control_restart(fp,c) 
controller c; 

FILE *fp; 

/************************************************/ 

{ complex_matrix ghat_tmp, phat_tmp, u_tmp, L, D; 
int i, k; 

ghat_tmp = cm_get2(fp); 
phat_tmp = cxn_get2 (fp); 

if ((ghat_tmp->nrows == c->ny) && (ghat_tmp->ncols == 1)) 

{ cm_assign(c->Ghat,ghat_tmp); 

> 

else 

{ printf("control_restart: Ghat isn't dimensioned right!") 
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> 

if ((phat_tmp->nrows == c->ny) && (phat_tmp->ncols «== c->nu)) 
{ cm_assign(c->Phat,phat_tmp); 

} 

else 

{ printf("control_restart: Phat isn't dimensioned right!"); 

} 

D = cm_get2(fp); 

L = cm_get2(£p); 

for (i = 0; i < c->nu + 1; i++) 

{ LDL_D(c->Qinv,i) = cv(D,i)->real; 
for (k = 0;k<i; k++) 

{ c_assign(LDL_L(c->Qinv,i,k),cm(L,i, k)); 

} 

} 

fscanf(fp,"sig2 %£, r Hf\n",&(c->sig2),6(c->r)) ; 

cm_free(D) ; 
cm_free(L) ; 
cmfree(ghat_tmp); 
cm_free(phat_tmp); 

/* set states that are not read in */ 
c->arm = 4; 

) t* control_restart */ 

/************************************************/ 

/* Controller save state for restart */ 

/* fp is open for writing. */ 

/************************************************/ 
void control_save(fp,c) 
controller c; 

FILE *fp; 

/************************************************/ 

{ int i, k, n; 

complex_matrix L, D; 

cm_put(fp/"Ghat" , c->Ghat); 
cm_put(fp, n Phat",c->Phat); 

n = (c->Qinv)->n; 

L = cm_new(n,n,0.0,0.0) ; 

D = cm_new(l,n,0.0,0.0); 

for (i = 0; i < c->nu + 1; i++) 

{ cv(D # i)->real = LDL_D(c->Qinv,i); 
for (k=0;k<i; k++) 

{ c_assign(cm(L,i,k),LDL_L(c->Qinv,i,k)); 

) 

) 

cm_put(fp,"Qinv diagonal",D); 
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cm_put(fp,"Qinv lower triangular",L) ; 

fprintf (fp, "sig2 %e, r 96e\n",c->sig2,c->r); 

cm_£ree(L); 
cm_£ree(D); 

} /* control_save */ 

/************************************************/ 

/* control algorithm */ 

/♦a**********************************************/ 

void control_iterate(y,u,c) 
float *u, *y; 
controller c; 

/♦a**********************************************/ 

{ int i, k, ny, nu, plant_changed; 

float rho, lambda, new_sig2, newr, maxu, certainty_index; 

ny ■ c->ny; 
nu = c->nu; 
lambda = c->l«unbda; 

/******************************************/ 

/* Map the measurements onto a complex y */ 
/******************************************/ 

for (i = 0; i < ny; i++) 

{ c_assign_real(cv(c->ym,i),y[2*i],y[2*i+l]); 

} 


/***********************************♦***/ 
/* Compute error and uncertainty ratio */ 
/* e = ym - Phat*um - Ghat */ 

/* eta = [1; urn] */ 

/* Qeta = inv(Qinv)*eta */ 

/* rho = real(eta'*Qeta) +1 */ 


/***************************************/ 

cm_mul(c->e,c->Phat,c->um); 
cm_sub(c->e,c->ym,c->e); 
cm_sub(c->e,c->e,c->Ghat); 

c_assign_real(cv(c->eta,0),1.0,0.0); 
cm_assign_part(c->eta,1,0, c->um, 0,0, nu,l); 

LDL_inv(c->Qeta,c->Qinv,c->eta); 
rho = cv_dotp(c->eta,c->Qeta) + 1.0; 

/* Check whether the plant has changed. */ 

/* Update uncertainty estimates and the */ 

/* inverse of Q. */ 

/****************************************/ 

new_sig2 = cv_norm2(c->e) / (rho * ((float) ny)); 
plant_changed = (new_sig2 > 20.0 * c->sig2) && (c->arm == 0) 
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i£ (plant_changed) 

{ 

/* Qinv = Qsmall + eta*eta' */ 

print£("\nplant change detected..."); 

LDL_eye(c->Qinvold, TOLER_CTL); 

LDL_update(c->Qinv, 1.0,c->Qinvold, 1.0,c->eta); 
c->arm = ARM_COUNT(c->nu); 

C->r - INITIAL_PERT; 
certainty_index = 0.0; 

) 

else 

{ 

/* Qinv = Qinvold + eta*eta' */ 

LDL_assign(c->Qinvold,c->Qinv); 

LDL_update(c->Qinv, 1.0,c->Qinvold, 1.0,c->eta); 
if (c->ann > 0) (c->arm) —; 
certainty_index = 

2.0 * sqrt(new_sig2*c->sig2) / (new_sig2 + c->sig2) / rho; 

) 

certainty_index = MINIMUM(1.0,3.0*certainty_index); 
c->sig2 += lambda*<new_sig2 - c->sig2); 

/***************************************/ 


/* Update plant estimates. */ 
/* Qeta = inv(Qinv)*eta */ 
/* Qetag = Qeta(0), Qetap * Qeta(l:nu) */ 
/* Ghat += e * Qetag' */ 
/* Phat += e * Qetap' */ 


LDL_inv(c->Qeta,c->Qinv # c->eta); 
c_assign(cv(c->Qetag,0),cv(c->Qeta / 0)); 
cm_assign_part(c->Qetap, 0,0, c->Qeta,l,0, nu,1); 

cv_xyH(c->Ghat,c->e,c->Qetag); 
cv_xyH(c->Phat,c->e,c->Qetap); 

/*******************************************/ 

/* Adjust perturbation size and direction. */ 

/* ipert cycles from 0 to 2*nu: */ 

/* 0 to nu-1 -- jog u(ipert) positive */ 

/* nu to 2nu-l -- jog u(ipert-nu) negative */ 

/* 2*nu — no perturbation */ 

/*******************************************/ 

newr = cv_norm(c->e) / (MIN_PN0RM + cv_norm(c->Phat)); 
c->r = MINIMUM(newr,10.0*c->r); 
c->r = MINIMUM(c->r,c->maxr); 

c->ipert = (c->ipert+l) % (2*nu+l); 


/★★★A***********************************/ 

/* Compute final control. */ 

/* u = -pseudoinv(Phat)*Ghat + pert */ 
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cm_inv2(c->Phatinv,c->Phat); 
cm_nnil(c->um,c->Phatinv,c->Ghat); 
cm_scale(c->um,-1.0,c->um); 

maxu c c->maxu; 

cv_limiter(c->um,(c->maxu * certainty_index)); 

c->r b c->r * certainty_index + (1.0-certainty_index) * c->maxr; 
if (c->ipert < nu) 

{ cm(c->\an / c->iport / 0)->real += c->r; 

} 

else if (c->ipert < 2*nu) 

{ cm(c->uxn / c->ipert-nu / 0)->real -= c->r; 

} 

/* Now assign the complex control to u */ 

for (i=0; i < nu; i++) 

{ u[2*i] = cv(c->um,i)->real; 
u[2*i+l] s cv(c->mn, i) ->imag; 

} 

} /* end control_iterate */ 


/************************************************^ 

/* Controller display. */ 

/************************************************/ 
void control_display(name, c) 
controller c; 
char *neune; 

{ static char *tf_str[2]; 
tf_str[FALSE] b "false"; 
tf_str[TRUE] = "true"; 

printf("\n\n==> Data for hd meas, Hd act controller \"%s\" <b=«, 
c->ny, c->nu, name); 

printf("\nsig2 b % e , r = He, ipert b Hd\n",c->sig2, c->r, c->ipert) 
cv_display("ym",c->ym); 
cm_display("Phat",c->Phat); 
cv_display("Ghat",c->Ghat); 

LDL_display("Qinv",c->Qinv); 
cv_display("urn",c->um); 

} /* end control_display) */ 
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Complex Matrix Utilities (cmplx6.h) 


,******************************************/ 
/* Simple math stuff */ 


/a*****************************************/ 

#undef C_TOLER 

#define C_TOLER (l.Oe-8) 

#undef MAXIMUM 

#define MAXIMUM(a,b) ((a) >= (b) ? (a) : (b)) 

#undef MINIMUM 

#define MINIMUM(a,b) ((a) <= (b) ? (a) : (b)) 

#undef ABS 

#define ABS(a) ((a) >= 0.0 ? (a) : -(a)) 

#undef LIMITER 

#define LIMITER(x # lim) MINIMUM(lim,MAXIMUM(x,-lim)) 

/******************************************/ 

/* Complex type and math. */ 

/******************************************/ 

typedef struct 
{ float real; 

float imag; 

} ‘complex; 

complex c_new(n,Re, Im) 
int n; 
float Re; 
float Im; 

{ complex c; 
int i; 

c ■ (complex) malloc(n*sizeof(*c)); 

for (i = 0; i < n; i++) 

{ (c[i]).real - Re; 

(c [i] ) . imag = Im; 

> 

return c; 

> 

void c_free(c) 
complex c; 

{ free(c); 

) 

void c_assign(to,from) 
complex to, from; 

{ (to)->real = (from)->real; (to)->imag = (from)->imag 

) 

void c_conj(z,x) 
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complex.z, x; 

{ (z)->real = (x)->real; 

(z)->imag = -((x)->imag); 

) 

void c_assign_real(to,Re,Im) 
complex to; 

£loat Re, Im; 

{ (to)->real = (Re); 

(to)->imag * (Im); 

} 

#de£ine c_zero(a) c_assign_real(a,0.0,0.0) 

void c_add(c,a,b) 
complex c,a,b; 

{ (c)->real = (a)->real + (b)->real; 

(c)->imag - (a)->imag + (b)->imag; 

) 

void c_sub(c,a,b) 
complex c, a, b; 

{ (c)->real ■ (a)->real - (b)->real; 

(c)->imag = (a)->imag - (b)->imag; 

) 

void c_mul(c,a # b) 
complex c, a, b; 

{ (c)->real = (a)->real * (b)->real - (a)->imag * (b)->imag; 

(c)->imag = (a)->real * (b)->imag + (a)->imag * (b)->real; 

} 

#define c_mag2(c) \ 

((c)->real * (c)->real + (c)->imag * (c)->imag) 

#define c_mag(c) \ 
sqrt(c_mag2(c)) 

/* Re{x*y'> */ 

#de£ine c_RexyH(x,y) \ 

((x)->real * (y)->real + (x)->imag * (y)->imag) 

/* z += x * conj(y) */ 
void c_xyH(z,x,y) 
complex z, x, y; 

{ (z)->real += (x)->real * (y)->real + (x)->imag * (y)->imag; 

(z)->imag += (x)->imag * (y)->real - (x)->real * (y)->imag; 

} 

/* z += x * y */ 
void c_xy(z,x,y) 
complex z, x, y; 

{ (z)->real += (x)->real * (y)->real - (x)->imag * (y)->imag; 

(z)->imag += (x)->imag * (y)->real + (x)->real * (y)->imag; 

/* z -= x * conj(y) */ 
void c_mxyH(z,x,y) 
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complex z, x, y; 

{ (z)-/real -■ (x)->real * (y)->real + (x)->imag * (y)->imag; 

(z)->imag -■ (x)->imag * (y)->real - (x)->real * (y)->imag; 

) 

/* z -* x * y */ 
void c_mxy(Z/X # y) 
complex z, x, y; 

{ (z)->real -= (x)->real * (y)->real - (x)->imag * (y)->imag; 

(z)->imag -= (x)->imag * (y)->real + (x)->real * (y)->imag; 

} 

/* z += a * x * y */ 
void c_axy(z,a,x,y) 
complex z, x, y; 
float a; 

{ (z)->real += (a)*((x)->real * (y)->real - (x)->imag * (y)->imag) 

(z)->imag += (a)*((x)->imag * (y)->real + (x)->real * (y)->imag) 

} 

/* z += a * x * yH */ 
void c_axyH(z # a,x,y) 
complex z, x, y; 
float a; 

{ (z)->real += (a)*((x)->real * (y)->real + (x)->imag * (y)->imag) 

(z)->imag += (a) *((x)->imag * (y)->real - (x)->real * (y)->imag) 

} 

/* z += a * x */ 
void c_ax(z,a,x) 
complex z, x; 
float a; 

{ (z)->real += (a)*((x)->real); 

(z)->imag += (a)*((x)->imag); 

} 

/* z = a*x */ 
void c_scale(z,a,x) 
complex z, x; 
float a; 

{ (z)->real = (a)*((x)->real); 

(z)->imag = (a)*((x)->imag); 

> 

/* prod += x'yz */ 
void c_xHyz(prod,x,y,z) 
complex prod, x, y, z; 

{ (prod)->real += 

(x)->real * ((y)->real*(z)->real - (y)->imag*(z)->imag) 

+ (x)->imag * ((y)->imag*(z)->real + (y)->real*(z)->imag); 
(prod)->imag += 

(x)->real * ((y)->imag*(z)->real + (y)->real*(z)->imag) 

(x)->imag * ((y)->real*(z)->real - (y)->imag*(z)->imag); 

} 

void c_div(z, x, y) 
complex z, x, y; 

{ float ymag2; 
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ymag2 = c_mag2(y); 

z->real = (x->real * y->real + x->imag * y->imag)/ymag2; 
z->imag = (x->imag * y->real - x->real * y->imag)/ymag2; 

) 

void c_display(title, c) 
char *title; 
complex c; 

{ printf("\n%s = (%13.6e,%13.6e)\n",title,c->real,c->imag); 

) 

/a**************************************************************/ 
/* Complex matrix and vector utilities */ 

/***************************************************************/ 

typedef struct 
{ int nrows; 
int ncols; 
complex ptr; 

} *complex_matrix; 

/* Allocate a matrix */ 

complexmatrix cm_new(n,m,Re,Im) 
int n, m; 
float Re, Im; 

{ complex_matrix z; 

z = (complex_matrix) malloc(sizeof(*z)); 
z->nrows = n; 
z->ncols = m; 

z->ptr = c_new(n*m,Re,Im); 
return z; 

} 

/* Deallocate a matrix */ 

void cm_free(z) 
complex_matrix z; 

{ free(z->ptr); 
free(z); 

) 

/* Get pointer to element (i,k) */ 

#define cm(z,i,k) \ 

( ((z)->ptr) + (i) + (k)*((z)->nrows) ) 

/* Get pointer to element (i) */ 

#define cv(z,i) \ 

(((z)->ptr) + (i)) 

/* z = x */ 

void cm_assign(z, x) 

complex_matrix z, x; 

{ int i,k; 

for (k = 0; k < x->ncols; k++) 

{ for (i = 0; i < x->nrows; i++) 
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{ c_as8ign(cm(z, i,k) ,cin(x, i,k)); 

) " 

} 

} 

/* z(iz+i,kz+k) = x(ix+i,kx+k) for i = 0:(nr-l), k = 0:(nc-l) 
void cm_assigm_part(z,iz,kz, x,ix,kx, nr,nc) 
complex_matrix z, x; 
int iz, kz, ix, kx, nr, nc; 

{ int i,k; 

for (k = 0; k < nc; k++) 

{ for (i = 0; i < nr; i++) 

{ c_assign(cm(z,iz+i,kz+k) ,cm(x,ix+i,kx+k)); 

} 

} 

} 

/* z(:,k) = v(:,i) */ 
cm_assign_col(z,k,v,i) 
complex_xnatrix z, v; 
int k,i; 

{ int r; 

for (r = 0; r < z->nrows; r++) 

{ c__assign (cm(z, r, k), cm(v, r, i) ) ; 

) 

) 

/* new z = x */ 
complex_matrix cm_copy(x) 
complex_matrix x; 

{ complex_matrix z; 

z = cm_new(x->nrows,x->ncols, 0.0,0.0); 
cm_assign(z,x); 
return z; 

} 

/* new z = x(ib:ie,kb:ke) */ 
complex_matrix cm_copy_part(x,ib, ie,kb,ke) 
complex_matrix x; 
int ib, ie, kb, ke; 

{ complex_matrix z; 
int nrow, ncol; 
nrow = (ie-ib+1); 
ncol = (ke-kb+1); 
z = cm_new(nrow,ncol, 0.0,0.0); 
cm_assign_part(z,0,0,x,ib,kb,nrow,ncol); 
return z; 

} 

/* z = 0 */ 
void cm_zero(z) 
complex_xnatrix z; 

{ int i,k; 

for (k = 0; k < z->ncols; k++) 

{ for (i = 0; i < z->nrows; i++) 

{ c_assign_real(cm(z,i,k),0.0,0.0); 

} 

} 
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/*»■!*/ 
void cm_eye(z) 
complex_matrix z; 

{ int i, n; 
cm_zero (z); 

n = MINIMUM(z->nrows,z->ncols); 

for (i = 0; i < n; i++) 

{ c_assign_real(cm(z,i,i),1.0,0.0); 

} 


complex_matrix cm_get(fp) 

FILE *fp; 

{ int nr, nc, i, k, readok; 
complex_matrix a; 

fscanf(fp, "%d%d n ,&nr,&nc); 

if ((nr < 1) || (nc < 1)) 

{ printf("\ncm_get: nrows or ncols < 1!"); 
return ((complex_matrix) NULL); 

> 

a = cm_new(nr,nc,0.0,0.0); 

for (i = 0; i < a->nrows; i++) 

{ for (k = 0; k < a->ncols; k++) 

{ readok = fscanf (fp, n s&e%e", & (cxn(a, i,k)->real) ,&(cm(a, i,k) - 
>imag)) ; 

if (readok == EOF) 

{ printf ( "cm_get: error reading element (s&d,%d) ! ", i,k) ; 
return ((complex_matrix) NULL); 

) 

) 

} 

return a; 

) 

void fskipline(fp) 

FILE *fp; 

{ char c; 

c = getc(fp); 
while (c != '\n') 

{ c = getc(fp); 

} 

> 

complex_matrix cm_get2(fp) 

FILE *fp; 

{ int nr, nc, i, k, readok; 
complex_matrix a; 
char line[80]; 
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£skipline(£p); 

fscanf(fp,"%d by %d",&nr,&nc); 

fskipline(fp); 

if {(nr < 1) || (nc < 1)) 

{ print£("\ncm_get: nrows or ncols < 1!"); 
return ((complex_matrix) NULL); 

} 

a = cm_new(nr,nc,0.0,0.0)? 

for (i = 0; i < a->nrows; i++) 

{ fBkipline(fp); 

for (k s 0; k < a->ncols; k++) 

{ readok = fscanf(fp,"%e %e",&(cm(a, i,k)->real),&(cm(a,i,k) 
>ixnag)); 

if (readok Is 2) 

{ printf ("cm_get: error reading element (%d,9&d) ! ", i,k) ; 
cm_free(a); 

return ((complex_matrix) NULL); 

) 

} 

fskipline(fp); 

} 

return a; 

} 

void cm_put(fp,title,a) 

FILE *fp; 
char *title; 
complex_matrix a; 

{ int i, k, 1; 

fprintf(fp,"Complex Matrix %s\n %d by %d",title, a->nrows, a- 
>ncols); 

for (i = 0; i < a->nrows; i++) 

{ 

fprintf (fp, "\nROW Ssd\n",i); 

for (k s o, 1 s 0; k < a->ncols; k++, 1++) 

if (1 == 2) 

{ 1 = 0; 

fprintf(fp,"\n"); 

) 

fprintf (fp," s&13.6e %13.6e ", 

cm(a,i,k)->real, cm(a,i,k)->imag); 

) 

fprintf(fp,"\n"); 


void cm_display(title,a) 
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char *title; 
complex_matrix a; 

{ int i, k, 1; 

printf("\nComplex Matrix %s",title); 

for (i = 0; i < a->nrows; i++) 

{ 

printf("\nROW %d\n",i); 

for (k = 0, 1 = 0; k < a->ncols; k++ # 1++) 

{ 

if (1 == 2) 

{ 1 = 0; 

printf("\n"); 

) 

printf("(%13.6e,%13.6e) ", 

cm(a,i,k)->real, cm (a,i,k)->imag); 

} 

} 

printf("\n"); 


void cv_display(title, a) 
char *title; 
complex_matrix a; 

{ int k, 1; 

printf("\nComplex Vector %s\n",title); 

for (k = 0, 1 = 0; k < (a->ncols)*(a->nrows); k++, 1++) 
{ 

if (1 == 2) 

{ 1 = 0; 

printf("\n"); 

> 

printf (" (S&13.6e,H13.6e) ", 

cv(a,k)->real, cv(a,k)->imag); 

) 

printf("\n"); 


/* c = a + b */ 
void cm_add(c,a,b) 
complex_matrix c, a, b; 

{ int i, k; 

for (k = 0; k < a->ncols; k++) 

{ for (i = 0; i < a->nrows; i++) 

{ c_add(cm(c,i,k),cm(a,i,k),cm(b,i,k)); 
) 

) 

} 

/* c = a - b */ 
void cm_sub(c,a,b) 
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cosqpleX'jnatrix c, a, b; 

< int i, k; 

for (k = 0/ k < a->ncols; k++) 

{ for (i e 0; i < a->nrows; i++) 

{ c_Bub(cm(c,i,k) / cm(a / i/k) f cm(b f i,k)); 

) 

) 

) 

/* c ■ a*b, c not the same as a or b */ 
void cm_mul(c,a,b) 
complex_matrix c, a, b; 

{ int i, k, 1; 

for (k = 0; k < b->ncols; k++) 

{ for (i = 0; i < a->nrows; i++) 

{ c_sero(cm(c,i,k)); 

for (1 = 0; 1 < b->nrows; 1++) 

{ c_xy(cm(c,i#k),cm(a,i # 1),cm(b,l,k)); 

} 

) 

} 

} 

/* cv_dotp = Re{v'*u> */ 
float cv_dotp(u, v) 
complex_matrix u, v; 

{ int i, n; 
float dotp; 

n = v->nrows * v->ncols; 
dotp = 0.0; 

for (i = 0; i < n; i++) 

dotp += c_RexyH{cv(u # i),cv(v,i)); 

return dotp; 

} 

/* ||v||-2 */ 
float cv_norm2(v) 
complex_matrix v; 

{ int i; 

float norm2; 

norm2 = 0.0; 

for (i = 0; i < (v->nrows)*(v->ncols); i++) 
{ norm2 += c_mag2(cv(v,i)); 

} 

return norm2; 

} 

/* ||v|| */ 

#define cv_norm(v) sqrt(cv_norm2(v)) 

/* z = a*x scale by a real number */ 
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void cm_scale(z,a,x) 
complexjmatrix z; 
float a; 

complexjmatrix x; 

{ int i; 

for (i = 0; i < z->nrows*z->ncols; i++) 

{ c_scale(cv(z,i),a,cv(x,i)); 

} 

> 

/* z += a*x */ 
void cv_ax(z,a,x) 
complex_matrix z; 
float a; 

complex_matrix x; 

{ int i; 

for (i e 0; i < z->nrows; i++) 

{ c_ax(cv(z # i),a,cv(x,i)); 

> 

} 

/* z b conj(x) */ 
void cv_conj(z,x) 
complexmatrix z # x; 

{ int i; 

for (i = 0; i < z->nrows*z->ncols 

; i++) 

{ c_conj(cv(z,i),cv(x,i)); 

> 

> 

/* x b (x)H */ 
void cm_H(x) 
complex_matrix x; 

{ int i,k,nr,nc; 
complex tmp; 
complex_matrix xtmp; 

nr b x->nrows; 
nc = x->ncols; 

if (nr b= l | nc == 1) /* vector case, <iuickly */ 

{ cv_conj(x,x); 
x->nrows = nc; 
x->ncols = nr; 

} 

else 

{ xtmp b cm_copy(x); /* matrix case, make copy, etc.. 

x->nrows = nc; 
x->ncols = nr; 
for (i b 0; i < nr; i++) 

{ for (k b 0; k < nc; k++) 

{ c_conj(cm(x,k,i),cm(xtmp,i,k)); 

> 

> 

cm_free(xtmp); 

} 

> 
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/* limits real and imag. parts of x to a certain range */ 
float cv_limiter(x,maxv) 
complex_matrix x; 
float maxv; 

/****************/ 

{ int i; 

float xr, xi, max_level, absxr, absxi; 
max_level = 0.0; 

for (i & 0; i < (x->nrows * x->ncols); i++) 

{ xr = cv(x,i)->real; 
xr c LIMITER(xr,maxv); 
absxr = ABS(xr); 

max_level = MAXIMUM(max_level,absxr); 

xi = cv(x,i)->imag; 
xi = LIMITER(xi,maxv); 
absxi = ABS(xi); 

max_level = MAXIMUM(max_level,absxi); 
c_assign_real(cv(x,i),xr,xi); 

) 

return max_level; 

} 

/* normalizes non-zero columns of x to unit magnitude */ 
void cm_normal(x) 
complex_matrix x; 

{ int i, k; 

float col_norm; 

for (k = 0; k < x->ncols; k++) 

{ col_norm = 0.0; 

for (i = 0; i < x->nrows; i++) 

{ col_norm += c_mag2(cm(x,i,k)); 

} 

col_norm = sqrt(col_norm); 

if (col_norm > 0.0) 

{ for (i * 0; i < x->nrows; i++) 

{ cm(x,i,k)->real /= col_norm; 
cm(x,i,k)->imag /= col_norm; 

) 

> 

} 

) 

/* z += x*yH outer product */ 
void cv_xyH(z,x,y) 
complexjnatrix z,x,y; 

{ int i, k; 

for (i = 0; i < z->nrows; i++) 

{ for (k = 0; k < z->ncols; k++) 

{ c_xyH(cm(z,i,k),cv(x,i),cv(y,k)); 

) 
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} 

} 

/* z += x*yT outer product */ 
void cv_xyT(z,x,y) 
complex_matrix z,x,y; 

{ int i, k; 

for (i = 0; i < z->nrows; i++) 

{ for (k = 0; k < z->ncols; k++) 

{ c_xy(cm(z # i # k),cv(x,i),cv(y,k)); 

) 

} 

} 

/* z -e x*yH outer product */ 
void cv_mxyH(z # X/y) 
complex_matrix z,x,y; 

{ int i, k; 

for (i = 0; i < z->nrows; i++) 

{ for (k = 0; k < z->ncols; k++) 

{ c_mxyH(cm(z,i,k),cv(x,i) # cv(y,k)); 

} 

} 

} 

/* z -= x*yT outer product */ 
void cv_mxyT(z,x / y) 
complex_matrix z,x f y; 

{ int i, k; 

for (i = 0; i < z->nrows; i++) 

{ for (k = 0; k < z->ncols; k++) 

{ c_mxy(cro(z,i # k),cv(x,i) , cv(y,k)); 

} 

} 

> 

/* z += a*x*yT outer product */ 
void cv_axyT(z,a,x,y) 
complex_matrix z; 
float a; 

complex_matrix x,y; 

{ int i, k; 

for (i = 0; i < z->nrows; i++) 

{ for (k = 0; k < z->ncols; k++) 

{ c_axy(cm(z,i,k),a,cv(x,i),cv(y,k)); 

} 

} 

} 

/* z += a*x*yH outer product */ 
void cv_axyH(z,a,x,y) 
complex_matrix z,x,y; 
float a; 

{ int i, k; 
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for (i = 0; i < z->nrows; i++) 

{ for (k e 0; k < z->ncols; k++) 

{ c_axyH(cm(z,i,k) ,a,cv(x,i),cv(y,k)); 

) 

} 

} 

/* z +* xT*y inner product, scalar result */ 
void cv_xTy(z,x,y) 
complex z; 
complex_matrix x,y; 

{ int if 

for (i ss 0; i < x->nrows; i++) 

{ c_xy(z,cv(y,i),cv(x,i)); 

) 

) 

/* z += xH*y inner product, scalar result */ 
void cv_xHy(z,x,y) 
complex z; 
complex_matrix x,y; 

{ int i; 

for (i = 0; i < x->nrows; i++) 

{ c_xyH(z,cv(y,i),cv(x,i)); 

> 

} 

/* z -= xT*y inner product, scalar result */ 
void cv_mxTy(z,x,y) 
complex z; 

complex_matrix x, y; 

{ int i; 

for (i = 0; i < x->nrows; i++) 

{ c_mxy(z,cv(y,i),cv(x,i)); 

} 


/* z -= xH*y inner product, scalar result */ 
void cvjnxHy(z,x,y) 
complex z; 

complex_matrix x, y; 

{ int i; 

for (i = 0; i < x->nrows; i++) 

{ c_mxyH(z,cv(y,i),cv(x,i)); 

} 


/* z += a*xT*y inner product, scalar result */ 
void cv_axTy(z,a,x,y) 
complex z; 
float a; 

complex_matrix x,y; 

{ int i; 
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for (i = 0; i < x->nrows; i++) 

{ c_axy(z,a,cv(y,i),cv(x,i)); 

) 

} 

/* z += a*xH*y inner product, scalar result */ 
void cv_axHy(z,a,x,y) 
complex z; 
float a; 

complex_matrix x,y; 

{ int i; 

for (i = 0; i < x->nrows; i++) 

{ c_axyH(z,a,cv(y,i),cv(x,i)); 

} 

} 

/*************************************************/ 


/* Utilities for Householder transformations. */ 
/a************************************************/ 

/* Finds a vector v having the same number of */ 
/* rows as a->ncols such that */ 
/* a(i,:) - 2*a(i,:)*v*vH is zero after */ 
/* column p. v is unit magnitude unless a is */ 
/* already in this form, when it is zero. */ 


complex_matrix house_v_row(a,i,p) 
complexmatrix a; 
int i, p; 

{ int m, n; 

complex_matrix v, x; 
complex phasorO; 
float xnorm, magxO; 

n = a->ncols; 
m = n-p; 


v = cm_new(n,1,0.0,0.0); 
if (m <= 0) 

{ return v; 

} 

/* Pull out X = a(i,p:(n-1))' */ 

x = cm_copy_part(a,i,i,p,(n-1)); 
cm_H(x); 

/* Get norm of x and |x(0)|, phasorO = x(0)/|x(0)|. 

xnorm = cv_norm(x); 
magxO = c_mag(cv(x,0)); 

phasorO = c_new(l,1.0,0.0); 
if (magxO > 0.0) 

{ phasor0->real = cv(x,0)->real / magxO; 
phasor0->imag = cv(x,0)->imag / magxO; 
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} 

/* Construct v in x, then assign. */ 

c_ax(cv(x,0),-xnorm,phasor0); /* x(0) =x(0) - xnorm*phasorO */ 
cm_normal(x); 

cm_assign_part(v,p,0, x,0,0, m, 1); 

c_free(phasor0); 
cm_£ree(x); 
return v; 


/* multiply x on left by 1 - 2*v*v' */ 

/* xbar may not be the same vector as x */ 
void house_Qx(xbar,v,x) 
complex_matrix xbar, v, x; 

{ complex_matrix vH, xtmp; 

if (v->nrows 1= x->nrows | xbar->nrows 1= x->nrows 
| xbar->ncols != x->ncols) 

{ return; 

} 

vH = cmcopy(v); 
cm_H(vH); 

xtmp = cm_new{l,x->ncols,0.0,0.0); 
cm_mul(xtmp,vH,x); 
cm_scale(xtmp,2.0,xtmp); 

cm_mul(xbar,v,xtmp); 
cmsub(xbar,x,xbar); 

cm_free(vH); 
cm_free(xtmp); 


/* multiply x on right by I - 2*v*v' */ 

/* xbar may not be the same vector as x */ 
void house_xQ(xbar,x,v) 
complex_matrix xbar, v, x; 

{ complexjmatrix vH, xtmp; 

if (v->nrows != x->ncols | xbar->nrows != x->nrows 
| xbar->ncols != x->ncols) 

{ return; 

} 

vH = cm_copy(v); 
cm_H(vH); 

xtmp = cm_new(x->nrows,1,0.0,0.0); 
cm_mul(xtmp,x,v); 
cm_scale(xtmp,2.0,xtmp); 

cm_mul(xbar,xtmp,vH); 
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cm_sub(xbar,x,xbar); 

cm_free(vH); 
cm_free(xtmp); 

} 

/*************************************************/ 
/* Utilities for an LDL' complex representation */ 
/*************************************************/ 

typedef struct 
{ int n; 

complex Lptr; 
float *Dptr; 

} *LDL_decomp; 

#define LDL_L(R,i,k) \ 

( ((R)->Lptr) + (i) + (k)*((R)->n) ) 

#define LDL_D(R,i) \ 

{(R)->Dptr)[i] 

/* Allocate an LDL decomposition */ 

LDL_decomp LDL_new(n) 
int n; 

{ LDL_decomp R; 
int i,k; 

R = (LDL_decomp) malloc(sizeof(*R)); 

R->n = n; 

R->Lptr = c_new(n*n,0.0,0.0) ; 

R->Dptr = (float (*)) malloc(n * sizeof(float)); 

for (i = 0; i < n; i++) 

{ LDL_D(R,i) = 0.0; 

c_assign_real(LDL_L(R,i,i),1.0,0.0); 

} 

return R; 

} 

void LDL_free(R) 

LDL_decomp R; 

{ c_free(R->Lptr); 
free(R->Dptr); 
free(R); 


void LDL_assign(Rto,Rfrom) 
LDL_decomp Rto, Rfrom; 

{ int i, k, n; 

n = Rfrom->n; 

for (i = 0; i < n; i++) 

{ LDL_D(Rto,i) = LDL_D(Rfrom,i); 
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for (k = 0; k < i; k++) 

{ classign(LDL_L(Rto,i,k) ,LDL_L(Rfrom, i,k)); 

} 

> 

> 

LDL_decomp LDL_get(fp) 

FILE *£p; 

{ int n, i, k, readok; 

LDL_decomp R; 

fscanf(fp,"%d",fcn); 

if (n < 1) 

{ printf("\nLDL_gets n < 1!"); 
return ((LDL_decomp) NULL); 

) 

R s LDL_new(n); 

for (k = 0; k < n; k++) 

{ readok « fscanf(fp,"%e n ,&LDL_D(R,k)); 
if (readok *= EOF) 

{ printf("LDL_get: error reading element %d of D!",k); 
return ((LDL_decomp) NULL); 

) 

} 

for (i = 0; i < n; i++) 

{ for (k = 0; k < n; k++) 

{ readok = fscanf (fp, "9&e%e n , 

&(LDL_L(R,i,k)->real),&(LDL_L(R,i,k)->imag)); 

if (readok == EOF) 

{ printf("LDL_get: error reading element (%d,%d) of L!",i,k) 
return ((LDL_decomp) NULL); 

} 

) 

> 

return R; 

} 

void LDL_display(title,R) 
char *title; 

LDL_decomp R; 

{ int i, k, 1; 

printf("\nLDL Decomposition %e” ,title); 

printf("\nDiagonal D:\n"); 

for (k = 0, 1 = 0; k < R->n; k++, 1++) 

{ 

if (1 == 5) 

{ 1 = 0; 

printf("\n"); 

} 
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printf("%13.6e ",LDL_D(R,k)); 

} 

printf("\nStrictly Lower Triangular Part of Complex Matrix L: n ) 
for (i = 1; i < R->n; i++) 

{ 

printf("\nROW %d\n",i); 

for (k = 0, 1 = 0; k < i; k++, 1++) 

{ 

if (1 == 2) 

{ 1 = 0; 

printf("\n"); 

) 

printf("{%13.ee,%13.6e) ", 

LDL_L(R,i,k)->real, LDL_L(R,i,k)->imag); 

) 

) 

printf("\n"); 

} 

/**************************************/ 

/* LDL update: Rnew = a*R + b*(x*x') */ 
/**************************************/ 

void LDL_update(R,a,Ro,b,x) 

LDL_decomp R, Ro; 
float a, b; 
complex_matrix x; 

{ int n, i, k; 
float nu, *rho; 
complex_matrix y; 
complex S, ctmp; 

y = cm_copy(x); 
n = Ro->n; 

/* Find y such that x = Lo * y; */ 

for (i = 0; i < n; i++) 

{ 

for (k = 0; k < i; k++) 

{ 

/* y(i) = y(i) - Lo(i,k)*y(k); */ 

c_mxy(cv(y,i),LDL_L{Ro,i,k),cv(y,k)); 

} 

} 

/* Now find rho and D. */ 

rho = f_new(n,0.0); 
nu = 1.0; 

for (i = 0; i < n; i++) 

{ 

LDL_D(R,i) = a * LDL_D(Ro, i) + b * nu * c_xnag2 (cv(y, i)) ; 

LDL_D(R,i) = MAXIMUM(0.0,LDL_D(R,i)); 


165 





ANC Final Report 


if '(LDL_D(R, i) ■■ 0.0) 

{ rhoti] = 0.0; 

) 

else 

{ rho[i] = b * nu / LDL_D(R,i); 

nu = a * nu * LDL_D(Ro,i) / LDL_D(R,i); 

} 

} 

/* Now compute L */ 

S = c_new(l,0.0,0.0); 
ctmp b c_new(l,0.0,0.0); 

for (i b l; i < n; i++) 

{ 

/* L(i,i-1) b Lo(i,i-l) + y(i)*y(i-1)'*rho[i-1]; */ 

c_assign(LDL_L(R,i,i-l),LDL_L(Ro,i,i-1)); 
c_axyH(LDL_L(R,i,i-l),rho[i-l],cv(y,i),cv(y,i-l)); 

c_assism_real(S, 0.0, 0.0); 

for (k = i-2; k >= 0; k—) 

{ 

/* S = S + Lo(i,k+l)*y(k+1); */ 

c_xy(S,LDL_L(Ro,i,k+l),cv(y,k+l)); 

/* L(i,k) = Lo(i,k) + (S + y(i))*y(k)'*rho(k); */ 

c_a s sign(LDL_L(R,i,k),LDL_L(Ro,i,k)>; 
c_add(ctmp,S,cv(y,i)); 

c_axyH(LDL_L(R,i,k),rho[k],ctmp,cv(y,k)); 

} 

) 

/* Free allocated resources. */ 

c_free(S); 
c_free(ctmp); 
f_free(rho); 
cm_free(y); 

) /* end LDL_update */ 

/A***********************************/ 

/* LDL inversion: u = inv(LDL')*v; */ 

void LDL_inv(u, R, v) 
complex_matrix u, v; 

LDL_decomp R; 

{ int i, k, n; 
float maxu2; 

n = R->n; 

/* Back substitute to get u = inv(L)*v */ 
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for (i = 0; i < n; i++) 

{ 

c_assign(cv(u,i),cv(v,i)); 

for (ks0;k<i; k++) 

{ 

c_nacy(cv(u, i) # cv(u,k) ,LDL_L(R, i,k)) ; 

} 

} 

/* Use pseudo-inverse of d */ 
maxu2 = c_mag2(cv(u,0)); 
for (i = 1; i < n; i++) 

{ maxu2 = MAXIMUM(maxu2,c_mag2(cv(u,i))); 

> 

for (i = 0; i < n; i++) 

{ if (c_mag2(cv(u, i)) > C_TOLER*maxu2) 

{ c_scale(cv(u, i) #(1.0/LDL_D(R,i) ) # cv(u,i)) 

) 

} 

/* Back substitute to get u */ 

for (i = (n-1); i >= 0; i—) 

{ 

for (k = (n-1); k > i; k--) 

{ 

c_mxyH(cv(u,i),cv(u,k),LDL_L(R,k,i)); 

} 

) 


/* LDL Pseudo Inverse constructor. */ 

/* R = dP2*eye + P'*P */ 


void LDL_pseudo(R,P,dP2) 

LDL_decomp R; 
complex_matrix P; 
float dP2; 

{ int i, k, n; 

LDL_decomp Rtmp; 
complex_matrix ptmp; 

n = R->n; 

Rtmp = LDL_new(n); 

ptmp = cm_new(n,1,0.0,0.0); 

for (i = 0; i < n; i++) 

{ LDL_D(Rtmp,i) = dP2; 

} 

/* Now accumulate the rows of P */ 
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for (i » 0; i < n; i++) 

{ 

for (k = 0; k < n; k++) 

{ c_conj(cv(ptmp,k),cm(P,i,k)); 

) 

LDL_update(R,1.0,Rtmp,1.0,ptmp); 
LDL_assign(Rtmp,R); 

} 

cm_free(ptmp); 

LDL_free(Rtmp); 


/* ainv * aH * inv(a * aH) */ 
void cm_inv2(ainv,a) 
complex_matrix ainv, a; 

{ complex_matrix u, v, aHinv; 

LDL_decomp aaH, aaHtmp; 
int nr, nc, i; 

nr = a->nrows; 
nc ■ a->ncols; 
aHinv = cm_copy(a); 

aaH = LDL_new(nr); 
aaHtmp = LDL_new(nr); 

u b cm_new(nr,1,0.0,0.0); 
v «= cm_new(nr, 1,0.0,0.0); 

/* Create aaH = a * aH as an LDL decomposition. */ 

for (i = 0; i < nc; i++) 

{ cm_assign_col(u,0,a,i); 

LDL_assign(aaHtmp,aaH); 

LDL_update(aaH,1.0,aaHtmp,1.0,u); 

) 

/* Back substitute aHinv(;,i) = inv(a*aH) * a(:,i) */ 

for (i=0; i<nc; i++) 

{ cm_assign_col(u,0,a,i); 

LDL_inv(v,aaH,u); 
cm_assign_col(aHinv,i,v,0); 

> 

cm_H(aHinv); 
cm_assign(ainv,aHinv); 

LDL_free(aaH); 

LDL_free(aaHtmp); 
cm_free(u); 
cm_free(v); 
cm_free(aHinv); 
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void LDL_eye(R,dxag_val) 
LDL_decomp R; 
float diag_val; 

{ int i,k; 

for (i = 0; i < R->n; i++) 
{ LDL_D(R,i) = diag_val; 

for (k = 0; k < i; k++) 

{ c_zero(LDL_L(R,i,k)); 

} 

) 

} 
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Floatingpoint Utilities (floatutl.h) 

/* floating point vector utilities */ 

float *£_new(n,fill_value) 
int n; 

float fill_value; 

{ int i; 
float *a; 

a = (float *) malloc(n*sizeof(float)); 
for (i =0; i < n; i++) 
a[i] = fill_value; 
return a; 

} 

void f_free(p) 
float *p; 

{ free(p); 

} 

void £_display(title,p,n) 
char *title; 
float *p; 
int n; 

{ int i,1; 

printf("\nflost vector %s\n",title); 

for (1=0, i = 0; i < n; i++) 

{ 

printf("%13.6e ",p(i]); 

1 ++; 

if (1 == 5) 1 = 0; 

} 

> 
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A/D Interface (adjib.h, adjib.c) 

ad_lib.h: 


/*.—7.*/ 

/* Function prototypes for ad_lib_gen */ 

/*- . - . -*/ 


void ad_setup(start_channel,nchannels); 
void ad_convert(float_buf _ptr,nchannels); 

ad_lib.c: 


/* .*/ 

/* File to define routines for ad control. */ 
/* L. Davis 5/9/96 */ 
/*.*/ 


/* type definitions for various word sizes */ 

typedef unsigned char onebyte; 
typedef unsigned short int twobyte; 
typedef unsigned int fourbyte; 

/* base addresses, etc. for 601 */ 

#define AD_BASE OxFOFAOOOO 
#define AD_DPR_LEN OxFFFF 
#define AD_CMD_OFFSET OxFFFO 
#define AD_WAKE_OFFSET OxFFEE 
#define AD_ZERO_CHAN 192 

twobyte *hststa = (twobyte *) (AD_BASE + 0xFFF6); 
twobyte *adstat = (twobyte *) (AD_BASE + 0xFFF4); 
twobyte *execcmd = (twobyte *) (AD_BASE + OxFFFC); 

/* routine addresses for the 601 */ 


#define DFSBUF 
#define DFDBUF 
#define SLSBUF 
#define SLWTEM 
#define SLDPXF 
#define XLCABF 
#define FASNSC 
fdefine STSNSC 
#define STSNSN 


((fourbyte) 0xl2E) 
((fourbyte) 0x134) 
((fourbyte) 0xl3A) 
((fourbyte) 0xl4C) 
((fourbyte) 0xl5E) 
((fourbyte) 0x17 0) 
((fourbyte) 0xl7C) 
((fourbyte) 0x188) 
((fourbyte) OxlBE) 


/* define low level tasks */ 


/*.-....*/ 

/* Allow permission to read/write */ 

/* code is in assembler at bottom. */ 

/* sys_permit(start_addr,#bytes); */ 

/*. ... */ 

static void sys_permit(); 
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#asm: 

sys_permit 

move.1 a2,-(sp) *save address 

movea.l d0,a2 *base address to a2 (argument 1) 

mo ve.1 dl,dO *offset to highest value (arg 2) 

moveq.l #7,dl 

os9 F$Permit *allow writing to these addresses/reading 

movea.l (sp)+ # a2 *restore a2 

rts 
#endasm 

/*...*/ 

/* Clear the DPR area out. */ 

/* ..*/ 

ad_clear_dpr() 

{ onebyte *bptr; 


for (bptr = (onebyte *) AD_BASE; 

bptr < (onebyte *) (AD_BASE + AD_DPR_LEN); bptr++) 
{ *bptr = (onebyte) 0; 

} 


/*. . —.—..*/ 

/* wait for acknowledgement of cmd */ 

/*....*/ 

void ad_wait_ack() 

{ int i = 0; 

while (!(0x0001 & *hststa)) /* check bit zero */ 

{ i++; i++; i—; i—; /* burn a few ops */ 

> 


void ad_reset_ack() 

{ *hststa = (*hststa & OxFFFE); 

} 

/*....*/ 

/* wait for afbO or afbO error */ 

/*—..*/ 

void ad_wait_afbO() 

{ int i = 0; 

while ((0x0050 6 *hststa) == 0) /* check bit zero */ 

{ i++; i++; i--; i--; /* burn a few ops */ 

} 

if ((0x0010 & *hststa) != 0) 

{ printf( n ad_wait_afb0: error in afbO execution!"); 

) 


/*--- . - -- 

/* wait for afbl or afbl error 
/*... 
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void »dlwait_afbl() 

{ int i = 0; 

while ((OxOOAO & *hststa) — 0) /* check bit zero */ 

{ i++; i++; i—; i—; /* bum a few ops */ 

} 

if ((0x0020 & *hststa) 1= 0) 

{ printf("ad_wait_afbl: error in afbl execution!"); 

) 

} 

/*.—...*/ 

/* execute command */ 

/*...*/ 

void ad_execute() 

{ *hststa = (*hststa & OxFFFE); /* reset ack */ 

*execcmd * (twobyte) 0; /* any write causes execution */ 

) 

/* ... */ 

/* transfer string to DPR cmd */ 

/* return new address */ 

/*--—.---*/ 

onebyte * ad_cmds_move(bptr,str) 
onebyte * bptr; 
char *str; 

{ int i; 

for (i = strlen(str)-1; i >= 0; i—) 

{ *(—bptr) ■ (onebyte) str[i]; 

) 

return bptr; 


/* ... */ 

/* transfer fourbyte data to DPR cmd. */ 
/* return new buffer address. */ 
/* alignment may not be even. */ 
/* .-.—..*/ 


onebyte *ad_cmd4_move(bptr,data) 
onebyte *bptr; 
fourbyte data; 

{ 

*(—bptr) = (onebyte) (data & OxOOOOOOFF) 
data = data >> 8; 

*(—bptr) = (onebyte) (data & OxOOOOOOFF) 
data = data >> 8; 

*(—bptr) = (onebyte) (data & OxOOOOOOFF) 
data = data >> 8; 

*(--bptr) = (onebyte) (data & OxOOOOOOFF) 
return bptr; 

) 
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/* . : .*/ 

/* View' memory byte by byte */ 

/* ..*/ 


void view_mem(base,offset,nbytes) 
onebyte *base; 
int offset/ nbytes; 

{ onebyte *ptr; 

static char *ctrl_str[32]; 
char c; 

ctrl_str[0] = "<nul>"; 
ctrl_str[l] «s "<soh>"; 
ctrl_str[2] = "<stx>"; 
ctrl_str[3] = "<etx>"; 
ctrl_str[4] = "<eot>"; 
ctrl_str[5] = "<enq>"; 
ctrl_str[6] = "<ack>"; 
ctrl_str[7] = "<bel>"; 
ctrl_str[8] = "<bs> 
ctrl_str[9] = "<ht> 
ctrl_str[10] = "<lf> 
ctrl_str[llj = "<vt> 
ctrl_str[12] = "<ff> 
ctrl_str[13] = "<cr> 
ctrl_str[14] = "<so> 
ctrl_str[15] = "<si> 
ctrl_str[16] = "<dle>"; 
ctrl_str[17] = "<dcl>"; 
ctrl_str[18] » "<dc2>"; 
ctrl_str[19] b "<dc3>"; 
ctrl_str[20] = "<dc4>"; 
ctrl_str[21] = *<nak>"; 
ctrl_str[22] = "<syn>"; 
ctrl_str[23] = "<etb>"; 
ctrl_str[24] = "<can>"; 
ctrl_str[25] = "<em> 
ctrl_str[26] = "<sub>"; 
ctrl_str[27] = "<esc>"; 
ctrl_str[28] = "<fs> 
ctrl_str[29] = "<gs> 
ctrl_str[30] = "<rs> 
ctrl_str[31] = "<us> 

printf("\n"); 

printf("\n+- 

printf("\n| MEMORY IMAGE 

printf (" \n+-+-+- 

printf("\n| address | hex | char 

printf (" \n+-+-+- 

for (ptr = base+offset; ptr < (base+offset+nbytes); ptr++) 
{ printf("\n| %8.8x | %2.2x ",(fourbyte) ptr,*ptr); 

c = *ptr; 
if (c < 32) 

{ printf(■| %5s |",ctrl_str[(int) c]); 

} 

else if (c >= 32 & c <= 126) 


■ + ■•) 
I") 
■ + '•) 
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{ printf("| He |",c); 

} : 

else if (c == 127) 

{ printf("| del |"); 

) 

else 

{ printf("|no char|"); 

} 

} 

printf (" \n+-+-+-+"); 

printf("\n"); 

) 


/* display hststa status */ 

/*....*/ 

void view_hststa() 

{ int bits[16], i; 

twobyte status_word, one = Oxl; 

status_word * *hststa; 

for (i = 0; i < 16; i++) 

{ bits[i] ■ ((one << i) & status_word) 1 = 0; 

} 

printf( n \n"); 

printf ("\n+-hststa status high byte- 

- +"); 

printf("\n| fund | funcO | fund | funcO | cntrl I not I not I 
downl |"); ' 

printf("\n| act. | act. | def. | def. | S | used | used I 
done |"); 

printf("\n| Hx | Hx | Hx | Hx |", 

bits[15],bits[14],bits[13],bits[12]); 
printf(" Hx | Hx | Hx | Hx |", 

bits[11],bits[10], bits[9],bits[8]); 

printf ( " \n+----hststa status low byte---— 

.+"); 

printf("\n| fund | funcO | fund | funcO | bad I cmd I downl I 
host |"); 

printf ("\n | coxnpl | compl | error | error | cmd j error | error | 
ack. |"); 

printf ("\n | Hx | Hx | Hx | Hx |", 

bits[7],bits[6],bits[5],bits[4]); 
printf (" Hx \ Hx | Hx | Hx |", 

bits[3],bits[2],bits[1],bits[0]); 

printf(-- 

.I"); 


.—.-.-...*/ 

/* display adstat semantics */ 


void view_adstat() 

{ twobyte status_word; 
int bits[16], i; 


175 












ANC Final Report 


twobyte one = Oxl; 
status_word = *adstat; 
for (i = 0; i < 16; i++) 

{ bits[i] = ((one << i) & status_word) != 0; 

} 


printf("\n"); 

printf ("\n+-adstat high byte 


printf("\n| exec. | not | pacer | not |wait on| local | not 
|use dbl|"); 

printf("\n| act. | used | def. | used | host | buffer| used 

buffer|"); 

printf ("\n I %x | %x | %x | %x |", 

bits[15],bits[14],bits[13],bits[12]); 
printf (" %x | %x | %x | | ", 

bits[11],bits[10] # bits[9],bits[8]); 

printf ("\n+-----adstat low byte ------------- — 

.+"); 

printf("\n| scan | conv. | not | buf 2 | buf 1 | buf 2 | buf 1 
which |"); 

printf("\n| done | done | used | def. | def. | rdy | rdy 

buf. |"); 

printf("\n| | %x | %x \ %x |", 

bits[7],bits[6],bits[5],bits[4]); 
printf (" %x | %x | %x \ %x |", 

bits[3],bits[2],bits[1],bits[0]); 

printf ( "\n+-- 

-| ", • 


I 


/*==========================================================*/ 

/* HIGHER LEVEL ROUTINES */ 

/*==========================================================*/ 

/*..*/ 

/* wake up the executive */ 

/*---.—.*/ 

void ad_wake_up() 

{ onebyte *bptr; 

/* allow DPR access */ 


sys_permit ((fourbyte) AD_BASE, (fourbyte) AD_DPR_LEN) ; 
ad_clear_dpr(); 

/* write wakeup commands into the DPR */ 

bptr = (onebyte *) (AD_BASE + AD_WAKE_OFFSET); 
bptr a ad_cmds_move(bptr,"EXEC ON\r n ); 

/* execute the command and try again. */ 

tsleep(20); 
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*hststa = (twobyte) 0; 

*adstat = (twobyte) 0; 

ad_execute(); 
tsleep(100); 

/* write wakeup commands into the DPR */ 

bptr ■ (onebyte *) (AD_BASE + AD_WAKE_OFFSET); 
bptr = ad_cmds_move(bptr,"EXEC ON\r"); 

/* execute the command and await acknowledgement */ 

tsleep(20); 

"hststa = (twobyte) 0; 

*adstat = (twobyte) 0; 

ad_execute(); 
tsleep(l); 

ad_wait_ack(); 

/* check for errors */ 

if (*hststa 1= 0x0001) 

{ printf("\nad_wake_up: hststa = (hex), wakeup 
failed!","hststa); 
return; 

) 


if ("adstat != 0x8000) 

{ printf("\nad_wake_up: adstat = %x (hex), wakeup 
failed!","adstat); 
return; 

) 

} 

f* .. */ 

/* Execute a single command to ad */ 

/*... */ 


void ad_do_.it (cmd_str) 
char *cmd_str; 

{ onebyte "bptr; 

bptr = (onebyte *) (AD_BASE + AD_CMD_OFFSET); 
bptr = ad_cmds_move(bptr,cmd_str); 
bptr = ad_cmds_move(bptr,"ENDC"); 

tsleep(10); 
ad_execute(); 
tsleep(l); 
ad_wait_ack(); 


/*.- . - . */ 

/* Set up a/d for scans */ 


177 




ANC Final Report 


/*.-.-...*/ 

void ad_setup(start_channel,nchannels) 
int start_channel, nchannels; 

{ onebyte *bptr; 

/* Wake up the a/d converter */ 

ad_wake_up(); 

/* Define "application function block (afb)" 0 */ 

bptr = (onebyte *) (AD_BASE + AD_CMD_OFFSET); 
bptr = ad_cmds_move (bptr, "BEGO") ; 

/* Define a single buffer in the dual-port memory for data. */ 

/* Set base offset zero, channel count. */ 

bptr = ad_cmd4_move(bptr,DFSBUF); 

bptr = ad_cmd4_move(bptr,(fourbyte) 0); 

bptr = ad_cmd4_move(bptr,(fourbyte) nchannels); 

/* Set single buffer transfers. */ 

bptr = ad_cmd4_move(bptr,SLSBUF); 

/* Make a/d wait until the data ready flag is reset by */ 

/* the host before doing the next conversion. */ 

bptr >s ad_cmd4_move (bptr, SLWTEM) ; 

/* Transmit data directly to the dual port RAM instead of */ 

/* to local a/d memory */ 

bptr = ad_cmd4_move(bptr,SLDPXF); 

/* Finish definition of AFBO, let the a/d digest it. */ 
bptr = ad_cmds_move(bptr,"ENDF"); 

bptr = ad_cmd4_move(bptr,(fourbyte) 0); /* insurance policy */ 

tsleep(20); 
ad_execute(); 
tsleep(1); 
ad_wait_ack(); 

/* Go ahead and process these setup commands once. */ 

ad_do_it("SELO"); /* select afbO */ 

ad_do_it("STRS"); /* execute once */ 

ad_wait_afbO(); /* wait until complete */ 

/* Redefine AFBO */ 

bptr = (onebyte *) (AD_BASE + AD_CMD_OFFSET); 
bptr = ad_cmds_move(bptr,"BEGO"); 

/* Use "fast" single scans. */ 
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bptr = ad_cmd4_move(bptr,FASNSC); 

bptr = ad_cmd4_move(bptr, ((fourbyte) AD_ZERO_CHAN + start_channel)) ; 
bptr = ad_cmd4_move(bptr, 

((fourbyte) (AD_ZERO_CHAN + start_channel + nchannels - 1))); 

/* Finish afbO definition, let a/d digest. */ 

bptr = ad_cmds_move(bptr,"ENDF"); 
bptr = ad_cmd4_move(bptr,(fourbyte) 0); 

tsleep(20); 
adexecute(); 
tsleep(l); 
ad_wait_ack(); 

/* now execute this continually until a software a/d reset, etc. */ 

ad_do_it("SELO"); 
ad_do_it("CONT")j 

/* Wait for a/d to get the message. */ 
ad_wait_afbO(); 

) 


/*..------ - */ 

/* Start conversion, wait and retrieve */ 

/* bufptr -- pointer to array of floats to put voltages in. */ 
/* nchannels — number of voltages to collect. */ 

/* --- ---*/ 

void ad_convert(bufptr,nchannels) 
float *bufptr; 
int nchannels; 

{ short int *adptr, i; 


/* Reset the data ready bit of adstat to start the scan. */ 

/* Wait for it to reassert when the scan is done. */ 

*adstat «t= (OxFD); 

while ((*adstat & 0x02) == 0) 

{ i = 0; i++; i--; /* burn some time */ 

) 

/* fetch the data from the buffer, converting to float voltages */ 

for (i = 0, adptr = (short int *) AD_BASE; i < nchannels; i++) 

{ *(bufptr++) = *(adptr++) / 3276.8; 

) 
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D/A Interface (dajib.h, dajib.c) 


da_lib.h: 

/*================================:==========*/ 

/* Function prototypes for PSI's da drivers */ 

/*BBBBBBSBB=B=BBBBSB=SBBBSBBBBSBBBBBB=BBB==B*/ 

void daO_setup(); /* call once to set up a board */ 
void dal_setup(); 
void da2_setup(); 

/* voltages — (float *) pointer to voltages between -<10.0 & 10.0 
*/ 

/* start_chan — (int) to which channel to write the first voltage 
*/ 

/* nchan -- (int) number of successive channels to write 

*/ 

void daO_convert(voltages,start_chan,nchan); . 
void dal_convert(voltages,start_chan,nchan); 
void da2_convert(voltages,start_chan,nchan); 

da_lib.c: 

/^======================================*/ 

/* PSI's da driver file. */ 

/*======================================*/ 


/* . — --- */ 

/* set da addresses */ 

,* - ---— . ---*/ 


#define DA0_BASE OxFOOOOOOO 
#define DA1_BASE OxFOOOOlOO 
#define DA2_BASE 0xF0000200 


#define CMD_OFFSET 0 
#define CHAN_OFFSET 2 
#define DATAOFFSET 4 
#define UPDATE_OFFSET 6 


static 

short 

int 

*da0_cmd = 

(DA0_BASE + 

static 

short 

int 

*dal_cmd = 

(DA1_BASE + 

static 

short 

int 

*da2_cmd = 

(DA2_BASE + 

static 

short 

int 

*da0_chan = 

(DA0_BASE + 

static 

short 

int 

*dal_chan = 

(DA1_BASE + 

static 

short 

int 

*da2_chan b 

(DA2_BASE + 

static 

short 

int 

*da0 data b 

(DA0_BASE + 

static 

short 

int 

*dal_data = 

(DA1_BASE + 

static 

short 

int 

*da2_data = 

(DA2_BASE + 

static 

short 

int 

*da0_update 

= (DA0_BASE 

static 

short 

int 

*dal_update 

b (DA1_BASE 

static 

short 

int 

*da2_update 

b (DA2 BASE 


CMD_OFFSET); 
CMD_OFFSET); 
CMD_OFFSET); 


CHAN_OFFSET); 
CHAN_OFFSET); 
CHAN_OFFSET); 

DATA_OFFSET); 
DATA_OFFSET); 
DATA_OFFSET); 

+ UPDATE_OFFSET); 
+ UPDATE_OFFSET); 
+ UPDATE_OFFSET); 


/* Set mode to immediate conversion on write-through */ 
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/* with, no auto-incrementing of the address register */ 
/* must write something to the update register to cause */ 
/* a conversion. */ 

#define NODE 0x0000 

/* Keep things in range */ 

#define MAXDA 2048 

#define DASENS (((float) MAXDA) / 10.0) 

#define damax(a,b) (((a) < (b)) ? (b) s (a)) 

#define damin(a,b) (((a) > (b)) ? (b) s (a)) 

#define da_encode(volts) (MAXDA + \ 

(short int) (DASENS * damin(10.0,damax(-10.0,volts)))) 


/* .. .- . */ 

/* system memory access permission routine. */ 

/* . — . ..—.. . */ 


static void sys_permit(); 
#asm: 

sys_permit 
move.1 
movea.1 
move.1 
moveq.1 
os9 

movea.1 
rts 

#endasm 


/*----- */ 

/* set up the d/a by writing each with the mode */ 

/*. .... .. */ 

void da0_setup() 

{ sys_permit ( (unsigned int) DA0_BA.SE, OxF) ; 

*da0_cmd = MODE; 

} 


void dal_setup() 

{ sys_permit((unsigned int) DA1_BASE, OxF); 
*dal_cmd = MODE; 

) 

void da2_setup() 

{ sys_permit((unsigned int) DA2_BASE, OxF); 
*da2_cmd = MODE; 

} 


/* —...-.*/ 

/* Send out a buffer full to daO for conversion. */ 
/* bufptr -- pointer to buffer of float voltages to */ 
/* write in [-10.0,10.0] */ 
/*....-.*/ 


void da0_convert(bufptr, start_chan, nchan) 

float *bufptr; 

int start_chan, nchan; 


a2, - (sp) 
d0,a2 
dl,d0 
#7, dl 

F$Permit 

(sp)+,a2 
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{ int i; 

float volts; 
short int value; 

/* Convert to an in-range integer for writing to da, then write */ 

for (i a start_chan; i < start_chan + nchan; i++) 

{ volts = *bufptr++; 

value a da_encode(volts); 

*daO_chan a i; 

*daO_data = value; 

) 

*daO_update =1; /* write anything to convert */ 


/*.-.*/ 

/* Send out a buffer full to dal for conversion. */ 
/* bufptr — pointer to buffer of float voltages to */ 
/* write in [-10.0,10.0] */ 
/*.......*/ 


void dal_convert(bufptr, start_chan, nchan) 

float *bufptr; 

int start_chan, nchan; 

{ int i; 

float volts; 
short int value; 

/* Convert to an in-range integer for writing to da, then write */ 

for (i = start_chan; i < start_chan + nchan; i++) 

{ volts = *bufptr++; 

value = da_encode(volts); 

*dal_chan = i; 

*dal_data = value; 

} 

*dal_update =1; /* write anything to convert */ 


/*” ... — .. . */ 

/* Send out a buffer full to da2 for conversion. */ 
/* bufptr — pointer to buffer of float voltages to */ 
/* write in [-10.0,10.0] */ 
/*...*/ 


void da2_convert(bufptr, start_chan, nchan) 

float *bufptr; 

int start_chan, nchan; 

{ int i; 

float volts; 
short int value; 

/* Convert to an in-range integer for writing to da, then write */ 

for (i = start_chan; i < start_chan + nchan; i++) 

{ volts = *bufptr++; 

value = da_encode(volts); 

*da2_chan = i; 

*da2_data = value; 

> 
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*da2_update =1; /* write anything to convert */ 

} 
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Appendix B. Electrical Schematics 


Analog Demodulator Cards 
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Sync Signal to Sine/Cosine Converter Cards 
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